\chapter{Input for {\dftbp}}


{\dftbp} can read the Human-friendly Structured Data format (HSD).  If you are
not familiar with the HSD format, a detailed description is given in appendix
\ref{sec:hsd}. The input file for {\dftbp} must be named \verb|dftb_in.hsd| and
must be present in the working directory. After processing the input, {\dftbp}
creates a file of the parsed input called \verb|dftb_pin.hsd|. This contains the
user input as well as default values for unspecified options.  The file also
contains the version number of the current input parser.  You should always keep
this file, since if you want to exactly repeat your calculation with a later
version of \dftbp{}, it is recommended to use this file instead of the original
input. (You must of course rename \verb|dftb_pin.hsd| into \verb|dftb_in.hsd|.)
This guarantees that you will obtain the same results, even if the defaults for
some non specified options have been changed in the meanwhile (as
\verb|dftb_pin.hsd| stores the original choice).

The following sections list properties and options that can be set in
{\dftbp} input. The first column of each of the tables of options
specifies the name of a property. The second column indicates the type
of the expected value for that property.  The letters ``l'', ``i'',
``r'', ``s'', ``p'', ``m'' stand for logical, integer, real, string,
property list and method type, respectively. An optional prefixing
number specifies how often (if more than once) this type must occur.
An appended ``+'' indicates arbitrary occurrence greater than zero,
while ``*'' allows also for zero occurrence.  Alternative types are
separated by ``|''.  Parentheses serve only to delimit groups of
settings.

Sometimes a property is only interpreted on meeting some condition(s).  If this
is the case, the the third column gives details of the requirement(s). The
fourth column contains the default value for the property.  If no default value
is specified (``-''), the user is required to assign a value to that property.
The description of the properties immediately follows the table.  If there is
also a more detailed description available for a given keyword somewhere else,
the appropriate page number appears in the last column.

Some properties are allowed to carry a modifier to alter the provided
value (e.g. converting between units). The possible modifiers are
listed between brackets ([]) in the detailed description of the
property. If the modifier is a conversion factor for a physical unit,
only the unit type is indicated (length, energy, force, time, etc.). A
list of the allowed physical units can be found in
appendix~\ref{app:units}.

\section{Main input}


The input file for {\dftbp} (\verb|dftb_in.hsd|)
must contain the following property definitions:
\begin{ptableh}
  \kw{Geometry} & p|m &  & - & \pref{sec:dftbp.Geometry} \\
  \kw{Hamiltonian} & m &  & - & \pref{sec:dftbp.Hamiltonian} \\
\end{ptableh}

Additionally optional blocks of definitions may be present:
\begin{ptableh}
  \kw{Driver} & m &  & \cb & \pref{sec:dftbp.Driver} \\
  \kw{Options} & p & & \cb & \pref{sec:dftbp.Options} \\
  \kw{Analysis} & p & & \cb & \pref{sec:dftbp.Analysis} \\
  \kw{ExcitedState} & p & & \cb & \pref{sec:dftbp.ExcitedState} \\
  \kw{ElectronDynamics} & p & & \cb & \pref{sec:dftbp.ElectronDynamics} \\
  \kw{REKS} & p & & \cb & \pref{sec:dftbp.REKS} \\
  \kw{InputVersion} & s & & None & \\
  \kw{ParserOptions} & p & & \cb & \pref{sec:dftbp.ParserOptions} \\
  \kw{Parallel} & p & & \cb & \pref{sec:dftbp.Parallel} \\
\end{ptableh}

\begin{description}
\item[\is{Geometry}] Specifies the geometry for the system to be
  calculated.  See p.~\pref{sec:dftbp.Geometry}.
\item[\is{Hamiltonian}] Configures the Hamiltonian and its options. See
  p.~\pref{sec:dftbp.Hamiltonian}.
\item[\is{Driver}] Specifies a geometry driver for your system.  See
  p.~\pref{sec:dftbp.Driver}.
\item[\is{Options}]  Various global options for the run. See
  p.~\pref{sec:dftbp.Options}.
\item[\is{Analysis}]  Post-run analysis and properties options. See
  p.~\pref{sec:dftbp.Analysis}.
\item[\is{ExcitedState}] Calculations in excited state of the system.
  See p.~\pref{sec:dftbp.ExcitedState}.
\item[\is{ElectronDynamics}] Real time electronic dynamics of the system.  See
  p.~\pref{sec:dftbp.ElectronDynamics}.
\item[\is{REKS}] Calculations in ensemble DFT for a strongly correlated system.
  See p.~\pref{sec:dftbp.REKS}.
\item[\is{InputVersion}] Program version, which the input file was written for.
\item[\is{ParserOptions}] Various options affecting the parser only.
  See p.~\pref{sec:dftbp.ParserOptions}.
\item[\is{Parallel}] Options affecting the MPI-parallel execution. See
  p.~\pref{sec:dftbp.Parallel}.
\end{description}


\section{Geometry}
\label{sec:dftbp.Geometry}

The geometry can be specified either directly by passing the
appropriate list of properties or by using the \is{GenFormat\cb}
method.
For molecular input the \is{xyzFormat\cb} is supported as well.
For periodic input the \is{VaspFormat\cb} supports reading POSCAR and CONTCAR files.

\subsection{Explicit geometry specification}

If the geometry is being specified explicitly, the following
properties can be set:

\begin{ptable}
  \kw{Periodic} & l & & No &  \\
  \kw{LatticeVectors} & 9r  & Periodic = Yes & - & \\
  \kw{TypeNames} & s+ &  & - &  \\
  \kw{TypesAndCoordinates}  & (1i3r)+  &  & - & \\
\end{ptable}
\begin{description}
\item[\is{Periodic}] Specifies if the system is periodic in all 3
  dimensions or is to be treated as a cluster. If set to \is{Yes},
  property \iscb{LatticeVectors} must be also specified.
\item[\is{LatticeVectors}]\modif{\modtype{length}} The $x$, $y$ and
  $z$ components of the three lattice vectors if the system is
  periodic.
\item[\is{TypeNames}] List of strings with the names of the elements,
  which appear in your geometry.
\item[\is{TypesAndCoordinates}] \modif{relative|\modtype{length}} For
  every atom the index of its type in the \is{TypeNames} list and its
  coordinates. If for a periodic system (\is{Periodic = Yes}) the
  modifier \is{relative} is specified, the coordinates are interpreted
  in the coordinate system of the lattice vectors.
\end{description}

Example: Geometry of GaAs:
\begin{verbatim}
Geometry = {
  TypeNames = { "Ga" "As" }
  TypesAndCoordinates [Angstrom] = {
    1  0.000000     0.000000     0.000000
    2  1.356773     1.356773     1.356773
  }
  Periodic = Yes
  LatticeVectors [Angstrom] = {
     2.713546     2.713546     0.
     0.           2.713546     2.713546
     2.713546     0.           2.713546
  }
}
\end{verbatim}

\subsection{GenFormat\cb}
\label{sec:dftbp.GenFormat}

You can use the generic format to specify the geometry (see
appendix~\ref{app:gen}). The geometry specification for GaAs would be
the following:
\begin{verbatim}
Geometry = GenFormat {
  2  S
  Ga As
  1 1     0.000000     0.000000     0.000000
  2 2     1.356773     1.356773     1.356773
  0.000000     0.000000     0.000000
  2.713546     2.713546     0.
  0.           2.713546     2.713546
  2.713546     0.           2.713546
}
\end{verbatim}
It is also possible to include the gen-formatted geometry from a file:
\begin{verbatim}
Geometry = GenFormat {
  <<< "geometry.gen"
}
\end{verbatim}

\subsection{xyzFormat\{\}}
\label{sec:dftbp.xyzFormat}

You can also use the xyz format to specify molecular geometries.
The geometry specification for a caffeine molecule would be the following:
\begin{verbatim}
Geometry = xyzFormat {
24
caffeine C8H10N4O2
C          1.07317        0.04885       -0.07573
N          2.51365        0.01256       -0.07580
C          3.35199        1.09592       -0.07533
N          4.61898        0.73028       -0.07549
C          4.57907       -0.63144       -0.07531
C          3.30131       -1.10256       -0.07524
C          2.98068       -2.48687       -0.07377
O          1.82530       -2.90038       -0.07577
N          4.11440       -3.30433       -0.06936
C          5.45174       -2.85618       -0.07235
O          6.38934       -3.65965       -0.07232
N          5.66240       -1.47682       -0.07487
C          7.00947       -0.93648       -0.07524
C          3.92063       -4.74093       -0.06158
H          0.73398        1.08786       -0.07503
H          0.71239       -0.45698        0.82335
H          0.71240       -0.45580       -0.97549
H          2.99301        2.11762       -0.07478
H          7.76531       -1.72634       -0.07591
H          7.14864       -0.32182        0.81969
H          7.14802       -0.32076       -0.96953
H          2.86501       -5.02316       -0.05833
H          4.40233       -5.15920        0.82837
H          4.40017       -5.16929       -0.94780
}
\end{verbatim}
It is also possible to include the xyz-formatted geometry from a file:
\begin{verbatim}
Geometry = xyzFormat {
  <<< "geometry.xyz"
}
\end{verbatim}

\subsection{VaspFormat\{\}}
\label{sec:dftbp.VaspFormat}

You can also use Vasp's \href{https://www.vasp.at/wiki/index.php/POSCAR}{POSCAR} or
\href{https://www.vasp.at/wiki/index.php/CONTCAR}{CONTCAR} format to specify periodic geometries.
The geometry specification for a molecular ammonia crystal would be the following:
\begin{verbatim}
Geometry = VaspFormat {
 H  N
    1.00000000000000
    5.01336000000000    0.00000000000000    0.00000000000000
    0.00000000000000    5.01336000000000    0.00000000000000
    0.00000000000000    0.00000000000000    5.01336000000000
 12 4
Cartesian
    2.19855889440000    1.76390058240000    0.88014548160000
    1.76390058240000    0.88014548160000    2.19855889440000
    0.88014548160000    2.19855889440000    1.76390058240000
    4.84115108400000    1.61941554720000    4.93981400880000
    4.35630903840000    2.49981169680000    3.63248012160000
    3.51957925440000    1.15357413600000    4.08403345680000
    4.08403345680000    3.51957925440000    1.15357413600000
    4.93981400880000    4.84115108400000    1.61941554720000
    3.63248012160000    4.35630903840000    2.49981169680000
    2.49981169680000    3.63248012160000    4.35630903840000
    1.15357413600000    4.08403345680000    3.51957925440000
    1.61941554720000    4.93981400880000    4.84115108400000
    1.37461317840000    1.37461317840000    1.37461317840000
    3.99815460000000    1.99105592400000    4.46364507600000
    4.46364507600000    3.99815460000000    1.99105592400000
    1.99105592400000    4.46364507600000    3.99815460000000
}
\end{verbatim}
It is also possible to include the POSCAR/CONTCAR formatted geometry from a file:
\begin{verbatim}
Geometry = VaspFormat {
  <<< "POSCAR"
}
\end{verbatim}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Driver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Driver}
\label{sec:dftbp.Driver}

The driver is responsible for changing the geometry of the input
structure during the calculation.\\ \bigskip Currently the following
methods are available:
\begin{description}
\item[\iscb{}] Static calculation with the input geometry.
\item[\iscb{SteepestDescent}] Geometry optimisation by moving atoms
  along the acting forces. See p.~\pref{sec:dftbp.SteepestDescent}.
\item[\iscb{ConjugateGradient}] Geometry optimisation using the
  conjugate gradient algorithm. See p.~\pref{sec:dftbp.ConjugateGradient}.
\item[\iscb{gDIIS}] Geometry optimisation using the modified gDIIS
  method. See p.~\pref{sec:dftbp.gDIIS}.
\item[\iscb{LBFGS}] Geometry optimisation using the LBFGS
  algorithm. See p.~\pref{sec:dftbp.LBFGS}.
\item[\iscb{FIRE}] Geometry optimisation using the FIRE
  algorithm. See p.~\pref{sec:dftbp.FIRE}.
\item[\iscb{SecondDerivatives}] Calculation of the second derivatives of the
  energy (the \index{Hessian}Hessian). See p.~\pref{sec:dftbp.SecondDerivatives}.
\item[\iscb{VelocityVerlet}] Molecular dynamics with the velocity
  Verlet algorithm. See p.~\pref{sec:dftbp.VelocityVerlet}.
\item[\iscb{Socket}] Hands over control to an external program via a socket
  interface. See p.~\pref{sec:dftbp.Socket}.
\end{description}


\subsection{SteepestDescent\cb}
\label{sec:dftbp.SteepestDescent}

\begin{ptable}
  \kw{MovedAtoms} & (i|s)+ &  & 1:-1 & \\
  \kw{MaxForceComponent} & r &  & 1e-4 & \\
  \kw{MaxSteps} &i &  & 200 & \\
  \kw{StepSize} &r &  & 100.0 & \\
  \kw{OutputPrefix} &s &  & "geo\_end" & \\
  \kw{AppendGeometries} & l & & \is{No} & \\
  \kw{Constraints} & (1i3r)* & LatticeOpt = No & \cb            & \\
  \kw{LatticeOpt}        & l & Periodic = Yes  & \is{No}        & \\
  \kw{FixAngles}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{FixLengths}        & 3l & FixAngles = Yes & \is{No No No} & \\
  \kw{Isotropic}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{Pressure}          & r & Periodic = Yes, LatticeOpt = Yes & 0.0
  & \\
  \kw{MaxAtomStep}       & r & MovedAtoms $\neq$ None\cb~ & 0.2 & \\
  \kw{MaxLatticeStep}    & r & Periodic = Yes, LatticeOpt = Yes & 0.2 & \\
  \kw{ConvergentForcesOnly} & l & SCC = Yes & \is{Yes} & \\
\end{ptable}
\begin{description}
\item[\is{MovedAtoms}] Indices of the atoms which should be moved. The
  atoms can be specified as a mixture of a list of atoms, ranges of
  atoms and/or the species of atoms. Index ranges are specified as
  \verb|start:end| (without white space as one word!), which
  inclusively selects all atoms between \verb|start| and \verb|end|.
  \index{Atom list}\index{List of atoms}
\begin{verbatim}
  MovedAtoms = 1:6
  # equivalent to MovedAtoms = { 1 2 3 4 5 6 }
\end{verbatim}
  Negative indices can be used to count backwards from the last atom
  (-1 = last atom, -2 = penultimate atom, etc.):
\begin{verbatim}
  MovedAtoms = 1:-1   # Move all atoms including the last
\end{verbatim}
  Species names can be used to select all atoms belonging to a given
  species:
\begin{verbatim}
  MovedAtoms = Ga  # select all Ga atoms
\end{verbatim}
  Various specifiers can be combined together:
\begin{verbatim}
  # Move atoms 1, 2, 3, all Ga atoms, and the last two atoms.
  MovedAtoms = 1:3 Ga -2:-1
\end{verbatim}

\item[\is{MaxForceComponent}]\modif{\modtype{force}} Optimisation is
  stopped, if the force component with the maximal absolute value goes
  below this threshold.

\item[\is{MaxSteps}] Maximum number of steps after which the optimisation should
  stop (unless already stopped by achieving convergence). Setting this value as
  -1 runs a huge() number of iterations.

\item[\is{StepSize}]\modif{\modtype{time}} Step size ($\delta t$)
  along the forces. The displacement $\delta x_i$ along the
  $i^\mathrm{th}$ coordinate is given for each atom as $\delta x_i =
  \frac{f_i}{2m}\delta t^2$, where $f_i$ is the appropriate force
  component and $m$ is the mass of the atom.

\item[\is{OutputPrefix}] Prefix of the geometry files containing the
  final structure.

\item[\is{AppendGeometries}] If set to \is{Yes}, the geometry file in
  the XYZ-format will contain all the geometries obtained during the
  optimisation (instead of containing only the last geometry).

\item[\is{Constraints}] Specifies geometry constraints. For every
  constraint the serial number of the atom is expected followed by the
  $x$, $y$, $z$ components of a constraint vector. The specified atom
  is not allowed to move along the constraint vector. If two
  constraints are defined for the same atom, the atom will only by
  able to move normal to the the plane containing the two constraining
  vectors.

  Example:
  \invparskip
\begin{verbatim}
  Constraints = {
    # Atom one can only move along the z-axis
    1  1.0  0.0  0.0
    1  0.0  1.0  0.0
  }
\end{verbatim}

\item[\is{LatticeOpt}] Allow the lattice vectors to change during
  optimisation. \is{MovedAtoms} can be optionally used with lattice
  optimisation if the atomic coordinates are to be co-optimised with
  the lattice.\footnote{This is functional but not very efficient at
    the moment.}

\item[\is{FixAngles}] If optimising the lattice, allow only the
  lengths of lattice vectors to vary, not the angles between them. For
  example if your lattice is orthorhombic, this option will maintain
  that symmetry during optimisation.

\item[\is{FixLengths}] If optimising the lattice with \is{FixAngles} =
  \is{Yes}, allow only the lengths of the specified lattice vectors to
  vary.

  Example:
  \invparskip
\begin{verbatim}
  Driver = ConjugateGradient {
    LatticeOpt = Yes
    FixAngles = Yes # Fix angles between lattice vectors
    FixLengths = {Yes Yes No} # Allow only lat. vector 3 to change length
  }
\end{verbatim}

\item[\is{Isotropic}] If optimising the lattice, allow only uniform
  scaling of the unit cell. This option is incompatible with
  \is{FixAngles}.

\item[\is{Pressure}]\modif{\modtype{pressure}} If optimising the lattice, set
  the external pressure, leading to a Gibbs free energy of the form $G = E + PV
  - TS$ being printed as well (the included entropy term is only the
  contribution from the electrons, therefore this is not the full free energy).

\item[\is{MaxAtomStep}] Sets the maximum possible line search step size
  for atomic relaxation.

\item[\is{MaxLatticeStep}] Sets the maximum possible line search step
  size for lattice optimisation. For \is{FixAngles} or \is{Isotropic}
  calculations this is as a fraction of the lattice vectors or the
  volume respectively.

\item[\is{ConvergentForcesOnly}] If using an SCC calculation, this
  option controls whether the geometry optimisation will prematurely
  stop (= \is{Yes}) if the SCC cycle does not converge at any
  geometric step.

\end{description}


\subsection{ConjugateGradient\cb}
\label{sec:dftbp.ConjugateGradient}

\begin{ptable}
  \kw{MovedAtoms} & (i|s)+ &  & 1:-1 & \\
  \kw{MaxForceComponent} & r & & 1e-4 & \\
  \kw{MaxSteps}          & i & & 200 & \\
  \kw{OutputPrefix}      & s & & "geo\_end" & \\
  \kw{AppendGeometries}  & l & & \is{No} & \\
  \kw{Constraints}       & (1i3r)* & & \cb & \\
  \kw{LatticeOpt}        & l & Periodic = Yes & \is{No} & \\
  \kw{FixAngles}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{Isotropic}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{Pressure}          & r & Periodic = Yes & 0.0 & \\
  \kw{MaxAtomStep}       & r & MovedAtoms $\neq$ None\cb~ & 0.2 & \\
  \kw{MaxLatticeStep}    & r & Periodic = Yes, LatticeOpt = Yes & 0.2 & \\
  \kw{ConvergentForcesOnly} & l & SCC = Yes & \is{Yes} & \\
\end{ptable}

See previous subsection for the description of the properties.

\subsection{gDIIS\cb}
\label{sec:dftbp.gDIIS}

\begin{ptable}
  \kw{Alpha} & r & & 0.1 & \\
  \kw{Generations} &i & & 8 & \\
  \kw{MovedAtoms} & (i|s)+ &  & 1:-1 & \\
  \kw{MaxForceComponent} & r & & 1e-4 & \\
  \kw{MaxSteps}          & i & & 200 & \\
  \kw{OutputPrefix}      & s & & "geo\_end" & \\
  \kw{AppendGeometries}  & l & & \is{No} & \\
  \kw{Constraints}       & (1i3r)* & & \cb & \\
  \kw{LatticeOpt}        & l & Periodic = Yes & \is{No} & \\
  \kw{FixAngles}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{Isotropic}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{Pressure}          & r & Periodic = Yes & 0.0 & \\
  \kw{MaxLatticeStep}    & r & Periodic = Yes, LatticeOpt = Yes & 0.2 & \\
  \kw{ConvergentForcesOnly} & l & SCC = Yes & \is{Yes} & \\
\end{ptable}

Specific properties for this method are:
\begin{description}
\item[\is{Alpha}] Initial scaling parameter to prevent the iterative
  space becoming exhausted (this is dynamically adjusted during the run).
\item[\is{Generations}] Number of generations to consider for the mixing.
\end{description}
See previous subsection for the description of the other
properties.\footnote{This approach is distinct from
  section~\ref{sec:dftbp.DIIS}, but uses a related algorithm based on
  Ref.~\cite{kovalenko-JCC-20-928} and comments from P.R.Briddon.}

\htcbsubsection{LBFGS}
\label{sec:dftbp.LBFGS}

\begin{ptable}
  \kw{Memory}            & i & & 20 & \\
  \kw{LineSearch}        & l & & No & \\
  \kw{MovedAtoms} & (i|s)+ &  & 1:-1 & \\
  \kw{MaxForceComponent} & r & & 1e-4 & \\
  \kw{MaxSteps}          & i & & 200 & \\
  \kw{OutputPrefix}      & s & & "geo\_end" & \\
  \kw{AppendGeometries}  & l & & \is{No} & \\
  \kw{Constraints}       & (1i3r)* & & \cb & \\
  \kw{LatticeOpt}        & l & Periodic = Yes & \is{No} & \\
  \kw{FixAngles}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{Isotropic}         & l & Periodic = Yes, LatticeOpt = Yes & \is{No} & \\
  \kw{Pressure}          & r & Periodic = Yes & 0.0 & \\
  \kw{MaxLatticeStep} & r & \parbox{0.4\textwidth}{LatticeOpt = Yes, LineSearch
    = Yes or setMaxStep = Yes} & 0.2 & \\
  \kw{MaxAtomStep}       & r & LineSearch = Yes or setMaxStep
  = Yes & 0.2 & \\
  \kw{ConvergentForcesOnly} & l & SCC = Yes & \is{Yes} & \\
\end{ptable}

Specific properties for this method are:
\begin{description}
\item[\is{Memory}] Number of last steps which are saved and used to calculate
  the next step via the lBFGS algorithm. The literature recommends that
  \is{Memory} should between 3 and 20~\cite{NoceWrig06}.
\item[\is{LineSearch}] Should a line search be performed, instead of a
  quasi-Newton step along the lBFGS direction.
\item[\is{SetMaxStep}] Limit maximum changes per iteration if the quasi-Newton
  step along the lBFGS direction is used.
\end{description}

\htcbsubsection{FIRE}
\label{sec:dftbp.FIRE}

The {\it fast inertial relaxation engine} of Bitzek {\it et
  al.}~\cite{Bitzek-PRL-97-170201}, using the default values from their paper
and a Leap-frog integrator. This can be faster than conjugate gradient and
competitive with LBFGS in some cases, but you should verify if FIRE is stable
and efficient for your system before using it over other optimisers.

\begin{ptable}
  \kw{TimeStep} & r & & 1.0 & \\
\end{ptable}

Specific properties for this method are:
\begin{description}
\item[\is{TimeStep}]\modif{\modtype{time}} Controls the maximum step size used
  in the dynamics integrator.
\end{description}

\subsection{SecondDerivatives\cb}
\label{sec:dftbp.SecondDerivatives}

Calculates the second derivatives\index{Hessian} of the energy (currently only
using a numerical differentiation of the forces). The derivatives matrix is
written out for the $i$, $j$ and $k$ directions of atoms $1 \ldots n$ as
\begin{equation*}
  \frac{\partial^2 E}{\partial x_{i1} \partial x_{i1}} \frac{\partial^2
    E}{\partial x_{j1} \partial x_{i1}} \frac{\partial^2 E}{\partial x_{k1}
    \partial x_{i1}} \frac{\partial^2 E}{\partial x_{i2} \partial x_{i1}}
  \frac{\partial^2 E}{\partial x_{j2} \partial x_{i1}} \frac{\partial^2
    E}{\partial x_{k2} \partial x_{i1}} \ldots \frac{\partial^2 E}{\partial
    x_{kn} \partial x_{kn}}
\end{equation*}
into {\it hessian.out}

\textbf{Note}: for supercell calculations, the derivatives are obtained at the $\mathbf{q}=0$ point,
irrespective of the k-point sampling used.

\textbf{Important:} In order to get accurate results for the second derivatives
(and the resulting frequencies) you must set a smaller self-consistent
tolerance than the default value in the
\islcb{Hamiltonian}{sec:dftbp.Hamiltonian} section. We suggest \is{SCCTolerance
  = 1e-7} or better. A less accurate tolerance can yield nonphysical vibrational
frequencies.

\begin{ptable}
  \kw{Atoms} & i+|m &  & 1:-1 & \\
  \kw{Delta} & r & & 1e-4 & \\
\end{ptable}

\begin{description}
\item[\is{Atoms}] Index of the atoms for which to calculate the second
  derivatives. The atoms can be specified via indices, index ranges
  and species. (See \is{MovedAtoms} in section \ref{sec:dftbp.SteepestDescent}.)
\item[\is{Delta}] Step size for numerical differentiation of forces to
  get the second derivatives of the energy with respect to atomic
  coordinates.
\end{description}


\subsection{VelocityVerlet\cb}
\label{sec:dftbp.VelocityVerlet}

The code propagates atomic motion using velocity Verlet dynamics with optional
thermostats or barostats to control the temperature and/or pressure. Information
is printed out during the simulation every
\kwl{MDRestartFrequency}{kw:dftbp.MDRestartFrequency} steps, and logged in the
file \verb|md.out| (see appendix \ref{sec:md.out}).

\begin{ptable}
  \kw{MovedAtoms} & (i|s)+ &  & 1:-1 & \\
  \kw{Steps} &i &  & - & \\
  \kw{TimeStep} & r & & - & \\
  \kw{KeepStationary} & l & & Yes & \\
  \kw{Thermostat} & m &  & - & \pref{sec:dftbp.Thermostat} \\
  \kw{OutputPrefix} &s & & "geo\_end" & \\
  \kw{MDRestartFrequency} & i & & 1 & \\
  \kw{Velocities} & (3r)* & & - & \\
  \kw{Barostat} & m & Periodic = Yes & - & \pref{sec:dftbp.Barostat}\\
  \kw{ConvergentForcesOnly} & l & SCC = Yes & \is{Yes} & \\
  \kw{Xlbomd} & p & XlbomdFast \textrm{not set} & & \pref{sec:dftbp.xlbomd} \\
  \kw{XlbomdFast} & p & Xlbomd \textrm{not set} & & \pref{sec:dftbp.xlbomd} \\
  \kw{Masses} & p & & & \pref{sec:dftbp.Masses} \\
  \kw{Plumed} & l & No & & \\
\end{ptable}

\begin{description}
\item[\is{MovedAtoms}] List of atoms to move during the MD. (See more
  detailed description on page \pref{sec:dftbp.SteepestDescent}.)

\item[\is{Steps}] Number of MD steps to perform. In the case of a
  thermostat using a \kwcb{TemperatureProfile} the number of steps is
  calculated from the profile.

\item[\is{KeepStationary}] Remove translational motion from the system.


\item[\is{TimeStep}]\modif{\modtype{time}} Time interval between two
  MD steps.

\item[\is{Thermostat}] Thermostating method for the MD simulation. See
  p.~\pref{sec:dftbp.Thermostat}.

\item[\is{OutputPrefix}] Prefix of the geometry files containing the
  final structure.

\item[\is{MDRestartFrequency}] Interval that the current geometry and
  velocities are written to the XYZ format geometry file and md.out
  (see section~\ref{sec:md.out}). In the case of \kw{SCC} MD runs, the
  charge restart information is also written at this interval
  overriding \is{RestartFrequency} (see
  section~\ref{kw:dftbp.RestartFrequency}).

\item[\is{Velocities}]\modif{\modtype{velocity}} Specified atomic velocities for
  all the atoms of the given structure (including ``velocities'' for any
  stationary atoms, which are silently ignored). This option can be used to
  restart an MD run, but make sure the geometry is consistent with the specified
  velocities. The easiest way to do this is to copy both from the same iteration
  of the XYZ file produced in the previous run (\textbf{Note:} the velocities
  printed in the XYZ files are specified in {\AA}/ps, so this should be set in
  the input). If restarting an \kw{SCC} MD run, we \textbf{strongly suggest} you
  use \is{ReadInitialCharges}, and in particular read charges for the geometry
  which you use to restart (iterations at which charges are written to disc are
  marked in the XYZ file, with the most recent being present in
  \verb|charges.bin| or \verb|charges.dat| depending on the option
  \is{WriteChargesAsText}).

\item[\is{Barostat}] Berendsen method barostat for the MD simulation. See
  p.~\pref{sec:dftbp.Barostat}.

\item[\is{ConvergentForcesOnly}] If using an SCC calculation, this option
  controls whether the molecular dynamics will prematurely stop (= \is{Yes}) if
  the SCC cycle does not converge at any geometric step. If the option is set to
  \is{False}, forces will be calculated using the non-converged charges and the
  molecular dynamics continues. In this case you should consider using
  \is{ForceEvaluation = 'Dynamics'} (or \is{ForceEvaluation = 'DynamicsT0'}) in
  the \is{DFTB} block as it gives more accurate forces for non-converged
  charges.

\item[\is{Xlbomd}] If present, extended Lagrangian type molecular dynamics
  is applied to speed up the simulation. For further options within the
  \is{Xlbomd} block see p.~\pref{sec:dftbp.xlbomd}.

\item[\is{Masses}] If present, over-ride the atomic masses from the Slater-Koster files. See
  p.~\pref{sec:dftbp.Masses}

\item[\is{Plumed}] Whether Plumed should be invoked in order to modify the
  forces and to drive a meta-dynamics.  The parameters of the meta-dynamics must
  be stored in the file \verb|plumed.dat| in the current directory. The file
  must be formatted according to Plumed's own format. (Note: This option
  requires a {\dftbp} binary built with Plumed support.)

\end{description}


\subsubsection{Thermostat}
\label{sec:dftbp.Thermostat}

\paragraph{\iscb{None}}

No thermostating during the run, only the initial velocities are set
according to either a given temperature or velocities, hence an
NVE\index{NVE ensemble} ensemble should be achieved for a reasonable
time step.
\begin{ptable}
  \kw{InitialTemperature} & r & & - & \\
\end{ptable}
\begin{description}
\item[\is{InitialTemperature}]\modif{\modtype{energy}} Starting
  velocities for the MD will be created according the
  Max\-well-Boltz\-mann distribution at the specified temperature.
  This is redundant in the case of specified initial velocities.
\end{description}

\paragraph{Andersen\cb}
\label{sec:dftbp.Andersen}

Andersen thermostat~\cite{andersen-JCP-72-2384} sampling an NVT\index{NVT ensemble} ensemble.

\textbf{Note:} Andersen thermostating has a reputation for suppressing diffusion and also prevents
accumulation of data for dynamical properties.

\begin{ptable}
  \kw{Temperature} & r|m & & - & \\
  \kw{ReselectProbability} &r & & - & \\
  \kw{ReselectIndividually} &l & & - & \\
  \kw{AdaptFillingTemp} & l & & No & \\
\end{ptable}
\begin{description}
\item[\is{Temperature}]\modif{\modtype{energy}} Target temperature of
  the thermostat.  It can be either a real value, specifying a
  constant temperature through the entire run or the
  \kwcb{TemperatureProfile} method specifying a changing temperature.
  (See p.~\pref{sec:dftbp.TemperatureProfile}.)
\item[\is{ReselectProbability}] Probability for re-selecting
  velocities from the Maxwell-Boltzmann distribution.
\item[\is{ReselectIndividually}] If \is{Yes}, each atomic velocity is
  re-selected individually with the specified probability. Otherwise
  all velocities are re-selected simultaneously with the specified
  probability.
\item[\is{AdaptFillingTemp}] If \is{Yes}, the temperature of the
  electron filling is always set to the current temperature of the
  thermostat. (The appropriate tag for the temperature of the electron
  filling is ignored.)

\end{description}

\paragraph{Berendsen\cb}
\label{sec:dftbp.Berendsen}

Berendsen thermostat~\cite{berendsen-JCP-81-3684} samples something
like an NVT\index{NVT ensemble} ensemble (but without the correct
canonical fluctuations, be aware of the ``flying ice cube'' problem
before using this thermostat~\cite{harvey-JCC-19-726}).

\begin{ptable}
  \kw{Temperature}      & r|m &                               & - & \\
  \kw{CouplingStrength} & r   & \is{Timescale} not set        & - & \\
  \kw{Timescale}        & r   & \is{CouplingStrength} not set & - & \\
  \kw{AdaptFillingTemp} & l   &                               & No & \\
\end{ptable}
\begin{description}
\item[\is{Temperature}]\modif{\modtype{energy}} Target temperature of
  the thermostat.  It can be either a real value specifying a constant
  temperature through the entire run or the \kwcb{TemperatureProfile}
  method specifying a changing temperature.  (See
  p.~\pref{sec:dftbp.TemperatureProfile}.)
\item[\is{CouplingStrength}] Dimensionless coupling strength for the
  thermostat (given as $\Delta t / \tau_t$ in the original paper where
  $\Delta t$ is the MD time step). Alternatively
  \is{Timescale}\modif{\modtype{time}} can be set directly as the
  characteristic length of time to damp the temperature towards the
  target temperature.  The \is{CouplingStrength} and \is{Timescale}
  options are mutually exclusive.

\item[\is{AdaptFillingTemp}] If \is{Yes}, the temperature of the
  electron filling is always set to the current temperature of the
  thermostat. (The appropriate tag for the temperature of the electron
  filling is ignored.)

\end{description}

\paragraph{NoseHoover\cb}
\label{sec:dftbp.NoseHoover}

Nos\'e-Hoover chain thermostat~\cite{martyna-mp-87-1117} sampling an
NVT\index{NVT ensemble} ensemble, currently with the chain coupled to
all of the atoms in the system.

\begin{ptable}
  \kw{Temperature} & r|m & & - & \\
  \kw{CouplingStrength} &r & & - & \\
  \kw{ChainLength} & i & & 3 & \\
  \kw{Order} & i & & 3 & \\
  \kw{IntegratorSteps} & i & & 1 & \\
  \kw{Restart} & m & & & \\
  \kw{AdaptFillingTemp} & l & & No & \\
\end{ptable}
\begin{description}
\item[\is{Temperature}]\modif{\modtype{energy}} Target temperature of
  the thermostat.  It can be either a real value, specifying a
  constant temperature through the entire run or the
  \kwcb{TemperatureProfile} method specifying a changing temperature.
  (See p.~\pref{sec:dftbp.TemperatureProfile}, but note that profiles are not
  well tested with this thermostat yet.)
\item[\is{CouplingStrength}]\modif{\modtype{Frequency}} Frequency of
  oscillation of the thermostating particles (see section~2.5 of
  Ref.~\cite{martyna-mp-87-1117}). This is typically related to the
  highest vibrational mode frequency of the system.
\item[\is{ChainLength}] Number of particles in the thermostat chain.
\item[\is{Order}] and \textbf{\is{IntegratorSteps}} See section~4.3 of
  Ref.~\cite{martyna-mp-87-1117}).
\item[\is{Restart}] Specifies the internal state of the thermostat
  chain, using three keywords (all three must be present): \kwcb{x},
  \kwcb{v} and \kwcb{g} containing the internal chain positions,
  velocities and forces respectively (these are provided in md.out).
  See also section~\ref{kw:dftbp.MDRestartFrequency}.
\item[\is{AdaptFillingTemp}] If \is{Yes}, the temperature of the
  electron filling is always set to the current temperature of the
  thermostat. (The appropriate tag for the temperature of the electron
  filling is ignored.)

\end{description}

\paragraph{TemperatureProfile\cb}
\label{sec:dftbp.TemperatureProfile}

Specifies a temperature profile during molecular dynamics. It takes as
argument one or more lines containing the type of the annealing
(string), its duration (integer), and the target temperature (real),
which should be achieved at the end of the given period. For example:
\begin{verbatim}
  Temperature [Kelvin] = TemperatureProfile {   # Temperatures in K
    constant        1    10.0   # Setting T=10 K for the 0th MD-step
    linear        500   600.0   # Linearly rising T in 500 steps up to T=600 K
    constant     2000   600.0   # Constant T through 2000 steps
    exponential   500    10.0   # Exponential decreasing in 500 steps to T=10 K
  }
\end{verbatim}
The annealing method can be \is{constant}, \is{linear} or
\is{exponential}, with the duration of each stage of the anneal
specified in steps of the driver containing the thermostat. The
starting temperature for each annealing period is the final target
temperature of the previous period, with the last step of each stage
being at the target temperature. Since the initial stage in the
temperature profile has no previous step, the default starting
temperature is set to 0.  In order to start the calculation from a finite
temperature for the first annealing period, a constant profile
temperature stage with a duration of one (or more) step(s) should be
specified as the first step (see the example).  The temperatures of
the stages are specified in atomic units, unless the \is{Temperature}
keyword carries a modifier, as in the example above.

\subsubsection{Barostat}
\label{sec:dftbp.Barostat}

Berendsen barostat~\cite{berendsen-JCP-81-3684} samples something like an
NPH\index{NPH ensemble} ensemble for MD (but without the correct
fluctuations). Options are provided for either isotropic or cell shape changing
pressure control. This can also be used in tandem with a thermostat
(p.~\pref{sec:dftbp.Thermostat}) for the NPT\index{NPT ensemble} ensemble. If the barostat
is used, a partial Gibbs free energy is reported in code output, of the form
\begin{equation*}
  G = E + PV - TS_\mathrm{electronic}.
\end{equation*}
This does not include structural entropy, only any electronic entropy. For
barostated constant energy simulations (no thermostat in use), the conserved
quantity is the sum of the kinetic and Gibbs partial energies.

\begin{ptable}
  \kw{Pressure}   & r & & - & \\
  \kw{Coupling}   & r & \is{Timescale} not set & - & \\
  \kw{Timescale}  & r & \is{Coupling} not set & - & \\
  \kw{Isotropic}  & l & & \is{Yes} & \\
\end{ptable}
\begin{description}
\item[\is{Pressure}]\modif{\modtype{pressure}} Sets the external
  target pressure.
\item[\is{Coupling}] Coupling strength for the barostat (given as
  $\beta \Delta t / \tau_p$ in the original paper where $\Delta t$ is
  the MD time step and $\beta$ the compressibility, so the coupling is
  technically dimensioned as reciprocal pressure, but this is usually
  ignored). Alternatively \is{Timescale}\modif{\modtype{time}} can be
  set directly ($\beta / \tau_p$) as the characteristic length of time
  to damp the pressure. Typically, $\beta$ is assumed to be either the
  experimental value or $\sim 1$, and
  Ref.~\cite{berendsen-JCP-81-3684} chooses the time scale to be
  around 10--100~fs. The \is{Coupling} and \is{Timescale} options are
  mutually exclusive.

\item[\is{Isotropic}] Should isotropic scaling of the unit cell be
  used, or can the cell shape vary? There is a slight inconsistency
  between the standard forms of these scalings in the literature,
  which is reproduced here, in brief the isotropic case scales the
  cell volume by a factor proportional to the differences in the
  instantaneous and expected pressures (i.e., the cube of the cell
  vectors), while the anisotropic case changes the cell vectors
  proportional to the difference instead.
\end{description}


\subsubsection{Extended Lagrangian Born-Oppenheimer dynamics}
\label{sec:dftbp.Extended Lagrangian Born-Oppenheimer dynamics}
\label{sec:dftbp.xlbomd}

For several systems Born-Oppenheimer molecular dynamics simulations can be
significantly sped up by using the extended Lagrangian formalism described in
Ref.~\cite{aradi-jctc-11-3357}. The XLBOMD integrator can be used in two
different modes:
\begin{itemize}
\item Conventional XLBOMD scheme (\is{Xlbomd}): The extended Lagrangian is used
  to predict the input charge distribution for the next time step, instead of
  taking charges that were converged for the geometry in the previous time
  step. The predicted starting charges should then require fewer SCC iterations
  to converge.

\item Fast XLBOMD scheme, \is{XlbomdFast} (one diagonalisation per time step):
  The extended Lagrangian is used to predict the population for each time
  step. This predicted population is then used to build the Hamiltonian, but in
  contrast to the conventional XLBOMD scheme, there is no self consistent cycle
  -- the forces are calculated immediately after the diagonalisation of the
  first Hamiltonian. The fast XLBOMD method usually only works for systems
  without SCC instabilities (e.g.\ wider gap insulators or molecules without
  degenerate states). See Ref.~\cite{aradi-jctc-11-3357} for details.
\end{itemize}

The extended Lagrangian dynamics can be activated by specifying either the
\is{Xlbomd} or the\\ \is{XlbomdFast} option block. Both methods offer the
following options:

\begin{ptable}
  \kw{IntegrationSteps} & i &  & 5 & \\
  \kw{PreSteps} & i &  & 0 & \\
\end{ptable}

\begin{description}
\item[\is{IntegrationSteps}] Number of time steps used for determining the
  population for the next time step. Currently, only integration schemes for 5,
  6 or 7 steps are implemented.
\item[\is{PreSteps}] Number of molecular dynamics time steps before the XLBOMD
  integration becomes activated.

  \textbf{Note:} At the step where the XLBOMD integrator becomes active, it is initialised with
  results of several subsequent converged MD steps, so a further \is{IntegrationSteps} + 1 steps
  will be carried out before the extended Lagrangian dynamics starts to predict the charges and
  accelerate the calculation.
\end{description}

The conventional \is{Xlbomd} block has the following specific options in
addition to the common XLBOMD settings above:

\begin{ptable}
  \kw{MinSccIterations} & i &  & 1 &  \\
  \kw{MaxSccIterations} & i & & 200 & \\
  \kw{SccTolerance} & r & & 1e-5 & \\
\end{ptable}

\begin{description}
\item[\is{MinSccIterations}] Minimum number of SCC iterations to perform at each
  time step during the extended Lagrangian dynamics.

\item[\is{MaxSccIterations}] Maximum number of SCC iterations to perform at each step in the
  extended Lagrangian dynamics. If this number of SCC iterations have been reached the forces will
  be calculated and the MD advances to the next time step. See the note in section~\ref{sec:dftbp.KLines}
  regarding non-convergent k-point sampling.

\item[\is{SccTolerance}] SCC convergence tolerance during the extended
  Lagrangian dynamics. Once this tolerance has been achieved the SCC cycle will
  stop and the forces will be calculated. You can use this parameter to override
  the \is{SccTolerance} parameter in the \is{DFTB} block for time steps where
  the extended Lagrangian integrator is active (This way, you can calculate
  populations with tight SCC tolerance when initialising the XLBOMD integrator,
  then use a less strict SCC tolerance once the integrator is active).
\end{description}

The \is{XlbomdFast} block has the following specific options in addition to
the common XLBOMD settings above:

\begin{ptable}
  \kw{TransientSteps} & i & & 10 & \\
  \kw{Scale} & r & & 1.0 & \\
\end{ptable}

\begin{description}
\item[\is{TransientSteps}] Enables a smoother transition between
  Born-Oppenheimer and extended Lagrangian dynamics by carrying out intermediate
  additional steps with full SCC convergence, during which the converged
  population and the one predicted by the extended Lagrangian integrator are
  averaged.

\item[\is{Scale}] Scaling factor for the predicted charge densities
  $\in(0,1]$. The optimal value is system dependent. One should take the highest
  possible value that still produces stable dynamics (good conservation of
  energy).
\end{description}

Example for conventional XLBOMD:
\begin{verbatim}
Xlbomd {
  IntegrationSteps = 6
  MinSccIterations = 2
  MaxSccIterations = 200
  SccTolerance = 1e-6
}
\end{verbatim}

Fast (SCC-free) XLBOMD with one diagonalisation per time step:
\begin{verbatim}
XlbomdFast {
  PreSteps = 5
  TransientSteps = 10
  IntegrationSteps = 5
  Scale = 0.5
}
\end{verbatim}

\textbf{Points to be aware of:}
\begin{itemize}
\item The extended Lagrangian (especially in the fast mode) needs special force
  evaluation giving more accurate forces for non-convergent charges. Therefore
  you must set the \is{ForceEvaluation} option to \is{'Dynamics'} (for
  simulations with finite electronic temperature) or to \is{'DynamicsT0'} (for
  simulations at 0~K electronic temperature) in the \is{DFTB} block (see
  p.~\pref{sec:dftbp.ForceEvaluation}).
\item The extended Lagrangian implementation only works for the $(N,V,E)$
  ensemble so far, so neither thermostats nor barostats are allowed.
\item The extended Lagrangian implementation currently cannot be used for
  spin-polarised or spin-orbit systems, or when electron filling methods other
  than \iscb{Fermi} filling (see p.~\pref{sec:dftbp.Fermi}) are used.
\end{itemize}

\subsubsection{Masses}
\label{sec:dftbp.Masses}

Provides values of atomic masses for specified atoms, ranges of atoms or chemical species. This is
useful for example to set isotopes for specific atoms in the system.

\begin{ptable}
  \kw{Mass} & p & & & \\
\end{ptable}

Any atoms not given specified masses will use the default values from the appropriate homonuclear
Slater-Koster file. An example is given below
\begin{verbatim}
  Masses {
    Mass {
      Atoms = 1:2
      MassPerAtom [amu] = 2.0
    }
  }
\end{verbatim}
where \kw{Atoms} specifies the atom or atoms which each have a mass of \kw{MassPerAtom} assigned.

\subsection{Socket\cb}
\label{sec:dftbp.Socket}

The code tries to connect to a socket interface to receive control instructions
from an external driver code.

\begin{ptable}
  \kw{File}  & s & \is{Host} not set & - & \\
  \kw{Prefix}& s & \is{Host} not set & ``/tmp/ipi\_'' for Protocol = \kwcb{i-PI}& \\
  \kw{Host}  & s & \is{File} not set & - & \\
  \kw{Port}  & i & \is{File} is set  & - & \\
  \kw{Verbosity} & i & & 0 & \\
  \kw{Protocol} & m & & \kwcb{i-PI} & \\
  \kw{MaxSteps} & i & & 200 & \\
\end{ptable}

\begin{description}
\item[\is{File}] Name of UNIX style file socket to connect to.
\item[\is{Prefix}] Prefix to the file name, in the case of i-PI dialogue, the
  defaults to the path and file start that i\_PI expects.
\item[\is{Host}] Name or ip address of internet host to connect to
  (``localhost'' also possible).
\item[\is{Port}] Port of host to connect to.
\item[\is{Verbosity}] Level of port traffic to document.
\item[\is{Protocol}] Choice of message protocol over the socket connection (only
  communication with \mbox{i-PI}~\cite{Ceriotti20141019} is currently
  supported).
\item[\is{MaxSteps}] Number of geometry steps before termination of the {\dftbp}
  instance.  Setting this value as -1 runs a huge() number of iterations.
\end{description}

Examples

First an ip address connection:
\invparskip
\begin{verbatim}
Driver = Socket {
    Host = localhost
    Port = 21012 # port number
    Verbosity = 0 # minimal verbosity
    Protocol = i-PI {} # i-PI interface
    MaxSteps = -1 # Run indefinitely
}
\end{verbatim}

Then a UNIX socket via the /tmp file system
\begin{verbatim}
Driver = Socket {
    File = "dftb" # The protocol defines a default path in this case
    Protocol = i-PI {} # i-PI interface
    MaxSteps = 1000 # Terminate this instance after 1000 steps
}
\end{verbatim}

Please note that this driver changes the default behaviour of some options to
remove (usually unneeded) file writing: \is{WriteDetailedOut} = \is{No} and
\is{WriteBandOut} = \is{No}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Hamiltonian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Hamiltonian}
\label{sec:dftbp.Hamiltonian}

Currently only a DFTB Hamiltonian is implemented, so you must set
\is{Hamiltonian = DFTB\cb}. The \iscb{DFTB} method may contain the
following properties:

\begin{ptable}
  \kw{SCC} &l &  & No  & \\
  \kw{SCCTolerance} &r& SCC = Yes & 1e-5 &  \\
  \kw{MaxSCCIterations} &i& SCC = Yes & 100 & \\
  \kw{EwaldParameter} & r & Periodic = Yes SCC = Yes & 0.0 & \\
  \kw{EwaldTolerance} & r & Periodic = Yes SCC = Yes & 1e-9 & \\
  \kw{ShellResolvedSCC} & l & SCC = Yes & No & \\
  \kw{Mixer} &m& SCC = Yes  & Broyden\{\} & \pref{sec:dftbp.Mixer} \\
  \kw{MaxAngularMomentum} &p&  & -  & \\
  \kw{Charge} &r& & 0.0 & \\
  \kw{SpinPolarisation} &m& SCC = Yes & \cb & \pref{sec:dftbp.SpinPolarisation} \\
  \kw{SpinConstants} &p& SpinPolarisation $\neq$ \cb~ & - &
  \pref{sec:dftbp.SpinConstants} \\
  \kw{ShellResolvedSpin} &l& ShellResolvedSCC = No & No & \\
  \kw{SpinOrbit} &m& SpinPolarisation $\neq$ Colinear\{\}~ & \cb &
  \pref{sec:dftbp.SpinOrbit} \\
  \kw{Solver} &m& & RelativelyRobust\cb  & \pref{sec:dftbp.Solver} \\
  \kw{Filling} &m&  & Fermi\cb & \pref{sec:dftbp.Filling} \\
  \kwcb{NonAufbau} & m & & - & \pref{sec:dftbp.NonAufbau} \\
  \kw{SlaterKosterFiles} &p|m&  & - & \pref{sec:dftbp.SlaterKosterFiles} \\
  \kw{OldSKInterpolation} & l & & No & \\
  \kw{PolynomialRepulsive} & p|m & & \cb & \\
  \kw{KPointsAndWeights} &(4r)+|m& Periodic = Yes & - &
  \pref{sec:dftbp.KPointsAndWeights} \\
  \kw{OrbitalPotential} &m & SpinPolarisation $\neq$ \cb~ & \cb &
  \pref{sec:dftbp.OrbitalPotential} \\
  \kw{ReadInitialCharges} &l& SCC = Yes & No & \\
  \kw{InitialCharges} & p & SCC = Yes & \cb & \\
  \kw{ElectricField} & p & SCC = Yes & \cb & \pref{sec:dftbp.ElectricField} \\
  \kw{Dispersion} & m & & \cb & \pref{sec:dftbp.Dispersion} \\
  \kw{HCorrection} & m & SCC = Yes & None \cb & \pref{sec:dftbp.hcorr}\\
  \kw{HalogenXCorr} & l & ThirdOrder(Full) = Yes, DftD3 & No & \pref{sec:dftbp.xcorr}\\
  \kw{ThirdOrder} & l & SCC = Yes & No & \\
  \kw{ThirdOrderFull} & l & SCC = Yes & No & \pref{sec:dftbp.DFTB3}\\
  \kw{RangeSeparated} & p &  & \is{None} & \pref{sec:dftbp.RangeSep} \\
  \kw{HubbardDerivs} & p & ThirdOrder(Full) = Yes  & - & \\
  \kw{OnSiteCorrection} & p & SCC = Yes & - & \pref{sec:dftbp.Onsites} \\
  \kw{Solvation} & m & & - & \pref{sec:dftbp.Solvation}\\
  \kw{Differentiation} &m&  & FiniteDiff & \pref{sec:dftbp.Differentiation} \\
  \kw{ForceEvaluation} &s& & "Legacy" & \\
  \kw{CustomisedHubbards} & p & SCC = Yes & & \\
  \kw{CustomisedOccupations} & p & SCC = Yes & & \\
  \kw{TruncateSKRange} & p & & & \\
  \hline
\end{ptable}

\begin{description}
\item[\is{SCC}] If set to \is{Yes}, a self consistent charge (SCC)
  calculation is made.

\item[\is{SCCTolerance}] Stopping criteria for the SCC.  Specifies the
  tolerance for the maximum difference in any charge between two SCC
  cycles.

\item[\is{MaxSCCIterations}] Maximal number of SCC cycles to reach
  convergence. If convergence is not reached after the specified
  number of steps, the program stops. However in cases where the
  calculation is not for a static structure (so \is{Driver} $\neq$
  \cb), this behaviour can be overridden using
  \kw{ConvergentForcesOnly}.

\item[\is{EwaldParameter}] Sets the dimensionless parameter $\alpha$
  in the Ewald electrostatic summation for periodic calculations. This
  controls the fraction of the Ewald summation occurring in real and
  reciprocal space. Setting it to zero or negative activates an
  automatic determination of this parameter (default and recommended
  behaviour). Setting it positive forces the code to use the supplied
  value. This is useful for very asymmetrical unit cells
  (typically a slab or nanowire with a huge surrounding vacuum region)
  since the memory demand of \dftbp{} can increase dramatically in
  these cases (due to storage of a long range real space neighbour
  list). To determine a suitable value of $\alpha$ for such a cell,
  you should initially reduce the vacuum region and run a test
  calculation, looking for the value of the automatically chosen Ewald
  parameter in the standard output. This is then a suitable choice for
  the original cell with the large vacuum region.

\item[\is{EwaldTolerance}] Sets the tolerance for Ewald summation in periodic
  systems.

\item[\is{ShellResolvedSCC}] If set to \is{Yes}, all distinct
  Hubbard $U$ values for the different atomic angular momenta shells
  are used, when calculating the SCC contributions. Otherwise, the
  value supplied for the $s$-shell is used for all angular
  momenta. Please note, that the old standard DFTB code was \emph{not}
  orbitally resolved, so that only the Hubbard $U$ for the $s$-shell
  was used. Please check the documentation of the SK-files you intend
  to use as to whether they are compatible with an orbitally resolved
  SCC calculation (many of the biological files do not use orbitally
  resolved charges), before you switch this option to \is{Yes}. Even
  if the Hubbard $U$ values for different shells are the same in the
  SK-files, this flag would still effect your results, since when it
  is set to \is{Yes}, any charge transfer between atomic shells will
  change the energy of the system compared to when it is set to set to
  \is{No}.

\item[\is{Mixer}]  Mixer type for mixing the charges in an SCC
  calculation. See p.~\pref{sec:dftbp.Mixer}.

\item[\is{MaxAngularMomentum}] Specifies the highest angular momentum
  for each atom type. All orbitals up to that angular momentum will be
  included in the calculation. Several main-block elements require
  $d$-orbitals, check the documentation of the SK-files you are using
  to determine if this is necessary. Possible values for the angular
  momenta are \is{s}, \is{p}, \is{d}, \is{f}.

  Example:
\begin{verbatim}
  MaxAngularMomentum = {
    Ga = "p"     # You can omit the quotes around the
    As = "p"     # orbital name, if you want.
  }
\end{verbatim}

  By using the \kw{SelectedShells} method it is also possible to
  combine shells from different Slater-Koster files together to treat
  atoms containing multiple shells with the same angular momentum. The
  shells to be picked from a particular atom type should be listed
  without any separating characters.  The list of shells of different
  atom types should be separated by white spaces.

  Example:
\begin{verbatim}
  # Defining sps* basis for Si and C by combining sp and s shells from
  # Si and Si2, and C and C2, resp.
  MaxAngularMomentum = {
    Si = SelectedShells { "sp" "s" }     # Si atom with sps* basis
    C = SelectedShells { "sp" "s" }      # C atom with sps* basis
  }

  # Note, that you have to modify the Slater-Koster file definition accordingly
  SlaterKosterFiles = {
    Si-Si = "Si-Si.skf" "Si-Si2.skf" "Si2-Si.skf" "Si2-Si2.skf"
    Si-C = "Si-C.skf" "Si-C2.skf" "Si2-C.skf" "Si2-C2.skf"
    C-Si = "C-Si.skf" "C-Si2.skf" "C2-Si.skf" "C2-Si2.skf"
    C-C = "C-C.skf" "C-C2.skf" "C2-C.skf" "C2-C2.skf"
  }
\end{verbatim}

  If for a given atomic type you pick orbitals from more than one species, you
  must specify an appropriate combinations of file names for the Slater-Koster
  tables in \islcb{SlaterKosterFiles}{sec:dftbp.SlaterKosterFiles}. For every
  atom type combination $n_{\text{SK1}}\times n_{\text{Sk2}}$ Slater-Koster
  files must be defined, where $n_{\text{SK1}}$ and $n_{\text{SK2}}$ are the
  number species combined to build up the shells of the two interacting
  atoms. The file names must be ordered with respect to the interacting species,
  so that the name for the second interacting species is changed first. Above
  you see an example, where an extended basis with an $s^*$-orbital was
  generated by introducing the new species "Si2" and "C2", containing the
  appropriate $s^*$-orbital for Si and C, resp., as only orbitals.

  In the case of multiple Slater-Koster files for a certain
  interaction, the repulsive data is read from the first specified
  file (e.g. \verb|Si-Si.skf|, \verb|Si-C.skf|, \verb|C-Si.skf| and
  \verb|C-C.skf| in the example above). The repulsive interactions in
  the other files are ignored. The mass for a certain species is read
  from the first SK-file for its homo-nuclear interaction.

  Non-minimal basis Slater-Koster data may be directly defined in the
  SK-files in future.

\item[\is{Charge}] Total charge of the system in units of the electron
  charge. Negative values mean an excess of electrons. If the keyword
  \is{FixedFermiLevel} is present (see section \ref{sec:dftbp.Filling}), then value
  specified here will be ignored.

\item[\is{SpinPolarisation}] Specifies if and how the system is spin
  polarised. See p.~\pref{sec:dftbp.SpinPolarisation}.

\item[\is{SpinConstants}] Specifies the atom type specific constants needed for
  the spin polarised calculations, in units of Hartree.  See
  p.~\pref{sec:dftbp.SpinConstants}.

\item[\is{SpinOrbit}] Specifies if the system includes Russel-Saunders
  coupling. See p.~\pref{sec:dftbp.SpinOrbit}

\item[\is{Solver}] Specifies which solver (eigensolver, Green's function, etc.)
  to use with the hamiltonian. See p.~\pref{sec:dftbp.Solver}.

\item[\is{Filling}] Method for occupying the one electron levels with
  electrons. See p.~\pref{sec:dftbp.Filling}.

\item[\is{SlaterKosterFiles}] Name of the Slater-Koster files for every atom
  type pair combination. See~\pref{sec:dftbp.SlaterKosterFiles}.

\item[\is{OldSKInterpolation}] If set to \is{Yes} (strongly
  discouraged), the look-up tables for the overlap and non-scc
  Hamiltonian contribution are interpolated with the same algorithm as
  in the \emph{old} {\dftb} code. Please note, that the new method
  uses a smoother function, is more systematic, and yields better
  derivatives than the old one.  This option is present only for
  compatibility purposes, and may be removed in the future.

\item[\is{PolynomialRepulsive}] Specifies for each interaction, if the
  polynomial repulsive function should be used. for every pairwise
  combination of atoms it should contain a logical value, where
  \is{Yes} stands for the use of a polynomial repulsive function and
  \is{No} for a spline. If a specific pair of species is not
  specified, the default value \is{No} is used.

  Example:
\begin{verbatim}
  # Use the polynomial repulsive function for Ga-Ga and As-As interactions
  # in GaAs
  PolynomialRepulsive = {
    Ga-Ga = Yes
    Ga-As = No
    # As-Ga unspecifed, therefore per default set to No
    As-As = Yes
  }
\end{verbatim}

  If you want to apply the same setting for all species pairs, you can
  specify the appropriate logical value as argument of the
  \is{SetForAll} keyword:

\begin{verbatim}
  # Using polynomial repulsive functions for all interactions in GaAs
  PolynomialRepulsive = SetForAll { Yes }
\end{verbatim}

\item[\is{KPointsAndWeights}]\modif{relative|absolute} Contains the
  special $k$-points to be used for the Bril\-louin-zone integration.
  See p.~\pref{sec:dftbp.KPointsAndWeights}. For automatically generated
  $k$-point grids the modifier should not be set.

\item[\is{OrbitalPotential}] Specifies which (if any) orbitally
  dependant contributions should be added to the DFTB energy and
  band structure. See p.~\pref{sec:dftbp.OrbitalPotential}.


\item[\is{ReadInitialCharges}] If set to \is{Yes} the first Hamiltonian is
  constructed by using the charge information read from the file
  \verb|charges.bin| or \verb|charges.dat| (depending on the option
  \is{WriteChargesAsText}, see section\ref{sec:dftbp.Options}).

\item[\is{InitialCharges}] Specifies initial charges, either for all atoms or
  for only selected ones. In order to specify it for all atoms, use the keyword
  \kw{AllAtomCharges} and supply the charge for every atom as a list of real
  values:
\begin{verbatim}
  InitialCharges = {
    AllAtomCharges = { # Specifies charge for each atom in an H2O molecule
       -0.88081627     # charge for atom 1 (O)
        0.44040813     # charge for atom 2 (H1)
        0.44040813     # charge for atom 3 (H2)
    }
  }
\end{verbatim}
  Alternatively you can specify charges individually on atoms or
  species using the \kw{AtomCharge} keyword. For every \is{AtomCharge}
  declaration you must provide a charge and the list of atoms,
  which should be initialised to that charge. (You can use the
  same format for the list of atoms, as described at the
  \is{MovedAtoms} keyword in the section for \isl{SteepestDescent}{sec:dftbp.SteepestDescent}):
\begin{verbatim}
  InitialCharges = { # Specifying charge for various species
    AtomCharge = {
      Atoms = H
      ChargePerAtom = 0.44040813
    }
    AtomCharge {
      Atoms = O
      ChargePerAtom = -0.88081627
    }
  }
\end{verbatim}
  Charge on atoms not appearing in any \is{AtomCharge}
  specification is set to be zero.


\item[\is{ElectricField}] Specifies an external electric field,
  arising currently from either an applied field or a distribution of
  electrostatic charges. See p.~\pref{sec:dftbp.ElectricField}.

\item[\is{Dispersion}] Specifies which kind of dispersion correction
  to apply. See p.~\pref{sec:dftbp.Dispersion}.

\item[\is{OnSiteCorrection}] Used to include the on-site matrix elements of
  Dom\'inguez {\it et al.}~\cite{dominguez15}. See p.~\pref{sec:dftbp.Onsites}.

\item[\is{Differentiation}] Specifies how to calculate finite difference
  derivatives in the force routines. See p.~\pref{sec:dftbp.Differentiation}.

\item[\is{ForceEvaluation}] Decides which expressions are used to calculate the ground state
  electronic forces.  See p.~\pref{sec:dftbp.ForceEvaluation}. \textbf{Note:} all methods give the same final
  forces when the charges are well converged. For non-converged charges however the resulting forces
  can differ considerably between methods.

\item[\is{CustomisedHubbards}] Enables overriding of the Hubbard U values for given species. If the
  option \is{OrbitalResolvedScc} has been set to \is{Yes}, you need to specify one Hubbard U value
  for each atomic shell in the basis of the given atom type, otherwise only one atomic value is
  required. For all species not specified in this block, the value(s) found in their respective
  Slater-Koster files will be used.

  \textbf{Warning:} This option is for experts only! Overriding values stored in the SK-files may
  result in \textbf{inconsistent results}. Make sure you understand the consequences when using this
  option.

  Example:
  \begin{verbatim}
  CustomisedHubbards {
    Si = 0.32 0.24
  }
  \end{verbatim}

\item[\kw{CustomisedOccupations}] Enables overriding the reference neutral atom
  electronic occupations of the shells. Note that the atom remains neutral since
  a corresponding ionic counter charge is implicitly added to its core. This
  option can be used, for example, to simulate effective doping by setting a
  slightly larger or smaller number of electrons on certain atoms.

  Example:
  \begin{verbatim}
  CustomisedOccupations{
    ReferenceOccupation{
      Atoms={1:30}
      p=2.01
    }
    ReferenceOccupation{
      Atoms={31:60}
      p=1.99
    }
  }
  \end{verbatim}

  The example above sets a filling population of +0.01e or -0.01e in the p shell
  of the corresponding atom indices. When the states are filled up, the electron
  excess or depletion results in a shift of the Fermi level in the bands.

  \textbf{Warning:} This option is for experts only! Overriding values stored in
  the SK-files may result in \textbf{inconsistent results}. Please look at the
  transport section of the
  \href{https://dftbplus-recipes.readthedocs.io/en/latest/}{dftb+ recipes} to
  see an example of the correct use of this option.

\item[\is{TruncateSKRange}] Enables overriding of the number of elements to be
  read from the Slater-Koster parameters, shortening the interaction range of
  atoms.

  \textbf{Warning:} This option is for experts only! Overriding values stored in
  the SK-files may result in \textbf{inconsistent results}. Make sure you
  understand the consequences when using this option.

  \begin{ptable}
    \kw{SKMaxDistance} & r &  & -- & \\
    \kw{HardCutOff} & l &  & No & \\
  \end{ptable}
  \begin{description}
    \item[\is{SKMaxDistance}]\modif{\modtype{length}} Length at which to cut the
      overlap and non-SCC interactions for all atoms in the system. If this
      length is longer than the distances in the Slater-Koster files it will
      have no effect.
    \item[\is{HardCutOff}] The Slater-Koster interpolation \dftbp{} produces
      will smoothly tail off to zero at a short distance beyond the table, which
      is controlled by \kw{OldSKInterpolation}. If \kw{HardCutOff} is set to
      \is{Yes}, then the distance set by \kw{SKMaxDistance} includes this tail,
      i.e., no interactions occur beyond that distance. In the case of \is{No}
      this zeroing tail extends slightly beyond the \kw{HardCutOff} distance.
  \end{description}

  Example:
  \begin{verbatim}
  TruncateSKRange = {
    SKMaxDistance [AA] = 4.0
    HardCutOff = Yes
  }
  \end{verbatim}

\end{description}


\subsection{Mixer}
\label{sec:dftbp.Mixer}

{\dftbp} currently offers the charge mixing methods \iscb{Broyden},
\iscb{Anderson}, \iscb{DIIS} (DIIS accel\-er\-ated simple mixer) and
\iscb{Simple} (simple mixer).

\subsubsection{Broyden\cb}
\label{sec:dftbp.Broyden}

Minimises the error function
\begin{equation*}
  E = \omega_0^2 \left| G^{(m+1)} - G^{(m)}\right| + \sum_{n=1}^m
  \omega_n^2 \left|
    \frac{n^{(n+1)} - n^{(n)}}{|F^{(n+1)}  - F^{(n)}|}
    + G^{(m+1)}
    \frac{F^{(n+1)}  - F^{(n)}}{|F^{(n+1)}  - F^{(n)}|} \right|^2
  \text,
\end{equation*}
where $G^{(m)}$ is the inverse Jacobian, $n^{(m)}$ and $F^{(m)}$ are
the charge and charge difference vector in iteration $m$. The weights
are given by $\omega_0$ and $\omega_m$, respectively. The latter is
calculated as
\begin{equation}
  \label{eq:omega_m}
  \omega_m = \frac{c}{\sqrt{F^{(m)}\cdot F^{(m)}}}
  \text,
\end{equation}
$c$ being a constant coefficient \cite{johnson-PRB-38-12807}.

The \iscb{Broyden} method can be configured using following properties:
\begin{ptable}
  \kw{MixingParameter} & r& & 0.2 & \\
  \kw{InverseJacobiWeight} & r & & 0.01 & \\
  \kw{MinimalWeight} & r & & 1.0 & \\
  \kw{MaximalWeight} & r & & 1e5 & \\
  \kw{WeightFactor} & r & & 1e-2 & \\
\end{ptable}
\begin{description}
\item[\is{MixingParameter}] Mixing parameter.
\item[\is{InverseJacobiWeight}] Weight for the difference of the
  inverse Jacobians ($\omega_0$).
\item[\is{MinimalWeight}] Minimal allowed value for the weighting
  factors $\omega_m$.
\item[\is{MaximalWeight}] Maximal allowed value for $\omega_m$.
\item[\is{WeightFactor}] Weighting factor $c$ for the calculation of
  the weighting factors $\omega_m$ in \eqref{eq:omega_m}.
\end{description}

Note: As the Broyden-mixer stores a copy of the mixed quantity for each SCC
iteration at a given geometry, you may consider to choose a different mixer with
lower memory requirements, if your system needs density matrix mixing (e.g.\
DFTB+U), is large and needs a high number of SCC-iterations
(\is{MaxSCCIteration}).


\subsubsection{Anderson\cb}
\label{sec:dftbp.Anderson}

Modified Anderson mixer \cite{eyert-JCP-124-271}.

\begin{ptable}
  \kw{MixingParameter} &r &  & 0.05 & \\
  \kw{Generations} &i &  & 4 & \\
  \kw{InitMixingParameter} & r & & 0.01 & \\
  \kw{DynMixingParameters} & (2r)* & & \cb  & \\
  \kw{DiagonalRescaling} & r & & 0.01 & \\
\end{ptable}
\begin{description}
\item[\is{MixingParameter}] Mixing parameter.
\item[\is{Generations}] Number of generations to consider for the
  mixing. Setting it too high can lead to linearly dependent sets of
  equation.
\item[\is{InitMixingParameter}] Simple mixing parameter used until the
  number of iterations is greater or equal to the number of
  generations.
\item[\is{DynMixingParameters}] Allows specification of different
  mixing parameters for different levels of convergence during the
  calculation. These are given as a list of tolerances and the mixing
  factor to be used below this level of convergence. If the loosest
  specified tolerance is reached, the appropriate mixing parameter
  supersedes that specified in \is{MixingParameter}.
\item[\is{DiagonalRescaling}] Used to increase the diagonal elements
  in the system of equations solved by the mixer. This can help to
  prevent linear dependencies occurring during the mixing
  process. Setting it to too large a value can prevent
  convergence. (This factor is defined in a slightly different way
  from Ref.~\cite{eyert-JCP-124-271}. See the source code for more
  details.)
\end{description}

Example:
\invparskip
\begin{verbatim}
  Mixer = Anderson {
    MixingParameter = 0.05
    Generations = 4
    # Now the over-ride the (previously hidden) default old settings
    InitMixingParameter = 0.01
    DynMixingParameters = {
      1.0e-2  0.1 # use 0.1 as mixing if more converged that 1.0e-2
      1.0e-3  0.3 # again, but 1.0e-3
      1.0e-4  0.5 # and the same
    }
    DiagonalRescaling = 0.01
  }
\end{verbatim}


\subsubsection{DIIS\cb}
\label{sec:dftbp.DIIS}

Direct inversion of the iterative space is a general method to
acceleration iterative sequences. The current implementation
accelerates the simple mix process.
\begin{ptable}
  \kw{InitMixingParameter} &r & & 0.2 & \\
  \kw{Generations} &i & & 6 & \\
  \kw{UseFromStart} &l & & \is{Yes} & \\
\end{ptable}
\begin{description}
\item[\is{MixingParameter}] Mixing parameter.
\item[\is{Generations}] Number of generations to consider for the mixing.
\item[\is{UseFromStart}] Specifies if DIIS mixing should be done right
  from the start, or only after the number of SCC-cycles is greater or
  equal to the number of generations.
\end{description}



\subsubsection{Simple\cb}
\label{sec:dftbp.Simple}

Constructs a linear combination of the current input and output
charges as $(1-x) q_{\text{in}}+ x q_{\text{out}}$.
\begin{ptable}
  MixingParameter & r & & 0.05 & \\
\end{ptable}
\begin{description}
\item[\is{MixingParameter}] Coefficient used in the linear
  combination.
\end{description}

\subsection{SpinPolarisation}
\label{sec:dftbp.SpinPolarisation}

In an SCC calculation, the code currently supports three different
choices for spin polarisation; non-SCC calculations are not spin
polarised.

\subsubsection{No spin polarisation (\is{\cb})}


No spin polarisation contributions to the energy or band-structure.

\subsubsection{Colinear\cb}
\label{sec:dftbp.Colinear}

Colinear spin polarisation in the $z$ direction.
%The initialisation of the calculation is spin restricted.
The following options can be
specified:
\begin{ptable}
  \kw{UnpairedElectrons} & r &  & 0  & \\
  \kw{RelaxTotalSpin} & l & & No & \\
  \kw{InitialSpins}      & p &  & \cb & \\
\end{ptable}
\begin{description}

\item[\is{UnpairedElectrons}] Number of unpaired electrons. This is kept
  constant during the run, unless the \is{RelaxTotalSpin} keywords says
  otherwise.

\item[\is{RelaxTotalSpin}] If set to \is{Yes}, a common Fermi-level is used for
  both spin channels, so that the total spin polarisation can change during
  run. In this case, the spin polarisation specified using the
  \is{UnpairedElectrons} keyword is only applied at initialisation. If set to
  \is{No} (default), the initial spin polarisation is kept constant during the
  entire run.

\item[\is{InitialSpins}] Optional initialisation for spin patterns. If
  this keyword is present, the default code behaviour is that the
  initial input charge distribution has a magnetisation of
  0. Otherwise if it is not present, the initial input charge
  distribution has a magnetisation matching the number of
  \is{UnpairedElectrons} keyword.

  The initial spin distribution for the input charges can be set by
  specifying the spin polarisation of atoms. For atoms without an
  explicit specification, a spin polarisation of zero is assumed. The
  \is{InitialSpins} property block must contain either the
  \kw{AllAtomSpins} keyword with a list of spin polarisation values
  for every atom, or one or more \kw{AtomSpin} blocks giving the spin
  for a specific group of atoms using the following keywords:
  \begin{ptable}
    \kw{Atoms} & (i|s)+ &  & -  & \\
    \kw{SpinPerAtom} & r &  & -  & \\
  \end{ptable}
  \begin{description}
  \item[\is{Atoms}] Atoms to specify an initial spin value. The atoms
    can be specified via indices, index ranges and species. (See
    \is{MovedAtoms} in section \ref{sec:dftbp.SteepestDescent}.)
  \item[\is{SpinPerAtom}] Initial spin polarisation for each atom in
  this \is{InitialSpins} block.
  \end{description}
  For atoms not appearing in any of the \is{SpinPerAtom} block, an
  initial spin polarisation of 0 is set.

Example (individual spin setting):
\invparskip
\begin{verbatim}
  SpinPolarisation = Colinear {
    UnpairedElectrons = 0.0
    InitialSpins = {
      AtomSpin = {
        Atoms = 1:2
        SpinPerAtom = -1.0
      }
      AtomSpin = {
        Atoms = 3:4
        SpinPerAtom = +1.0
      }
    }
  }
\end{verbatim}

Example (setting all spins together):
\invparskip
\begin{verbatim}
  SpinPolarisation = Colinear {
    UnpairedElectrons = 0.0
    InitialSpins = {
      AllAtomSpins = { -1.0 -1.0 1.0 1.0 } # Atoms 1,2: -1.0, atoms 3,4: 1.0
    }
  }
\end{verbatim}
\end{description}

\subsubsection{NonColinear\cb}
\label{sec:dftbp.NonColinear}

Non-collinear spin polarisation with arbitrary spin polarisation vector
on every atom. The only option allowed is to set the initial spin
distribution:
\begin{ptable}
  \kw{InitialSpins}       & p &  & \cb & \\
\end{ptable}
\begin{description}
\item[\is{InitialSpins}] Initialisation of the $x$, $y$ and $z$
  components of the initial spins for atoms. The default code
  behaviour is an initial spin polarisation of (0 0 0) for every atom.

  The initial spin distribution can be set by specifying the spin
  polarisation vector on all atoms using the \kw{AllAtomSpins}
  keyword or by using one or more \kw{AtomSpin} blocks with the
  following options:
  \begin{ptable}
    \kw{Atoms} & (i|s)+ &  & -  & \\
    \kw{SpinPerAtom} & (3r)+ &  & -  & \\
  \end{ptable}
  \begin{description}
  \item[\is{Atoms}] Atoms to specify an initial spin vector. The atoms
    can be specified via indices, index ranges and species. (See
    \is{MovedAtoms} in section \ref{sec:dftbp.SteepestDescent}.)
  \item[\is{SpinPerAtom}] Initial spin polarisation for each atom in
    this \is{InitialSpins} block.
  \end{description}
  For atoms not appearing in any of the \is{SpinPerAtom} block, the
  vector (0 0 0) is set.

  Please note, that in contrast to the collinear case, in the
  non-collinear case you must specify the spin vector (3 components!)
  for the atoms.

  Example:
  \invparskip
\begin{verbatim}
  SpinPolarisation = NonColinear {
    InitialSpins = {
      # Setting the spin for all atoms in the system
      AllAtomSpins = {
        0.408 -0.408  0.816
        0.408 -0.408  0.816
       -0.408  0.408 -0.816
       -0.408  0.408 -0.816
      }
    }
  }
\end{verbatim}

  Example:
\invparskip
\begin{verbatim}
  SpinPolarisation = NonColinear {
    InitialSpins = {
      AtomSpin = {
        Atoms = 1:2
        SpinPerAtom = 0.408 -0.408 0.816
      }
      AtomSpin = {
        Atoms = 3:4
        SpinPerAtom = -0.408  0.408 -0.816
      }
    }
  }
\end{verbatim}
\end{description}

\subsubsection{SpinConstants}
\label{sec:dftbp.SpinConstants}

This environment suplies the atomic constants required for either spin polarised
calculations or when evaluating properties which depend on spin interactions
(triplet excitations for example).  In these cases, for each atomic species in
the calculation the spin coupling constants for that atom must be specified.

When \is{ShellResolvedSCC} = \is{No}, an extra keyword in this block controls
whether the spin constants are resolved by shell or are identical for all
shells: \kw{ShellResolvedSpin}, defaulting to the same value as
\is{ShellResolvedSCC}.

When shell resolved spin constants are specified, they must be ordered with
respect to the pairs of shells they couple, such that the index for the second
shell increases faster. For an $spd$-basis, that gives the following ordering:
\begin{equation*}
  w_{ss}, w_{sp}, w_{sd}, \dots,
  w_{ps}, w_{pp}, w_{pd}, \dots,
  w_{ds}, w_{dp}, w_{dd}, \dots
\end{equation*}

Example (GGA parameters for H$_2$O):
\begin{verbatim}
  SpinConstants = {
    O = {
      # Wss  Wsp    Wps    Wpp
      -0.035 -0.030 -0.030 -0.028
    }
    H = {
      # Wss
      -0.072
    }
  }
\end{verbatim}

Several standard values of atomic spin constants are given in
appendix~\ref{app:spinconst}. Constants calculated with the same
density functional as the SK-files should be used. This input block
may be moved to the SK-data definition files in the future.

When using the \is{SelectedShells} method for the keyword
\is{MaxAngularMomentum}, the spin constants are listed as an array of
values running over $\text{SK1}\text{SK2}\ldots$ in the same order as
listed for SlaterKosterFiles.

\begin{verbatim}
  SpinConstants = { # not real values, only an example
    Si = {
      # Wss  Wsp    Wss*
      -0.035 -0.030 -0.01
      # Wps  Wpp    Wps*
      -0.030 -0.037 -0.02
      # Ws*s Ws*p   Ws*s*
      -0.01  -0.02  -0.01
    }
\end{verbatim}

For cases where \is{ShellResolvedSpin} = \is{No}, the spin constant for the the
highest occupied orbital of each atom should be supplied: Example (GGA parameters
for H$_2$O):
\begin{verbatim}
  SpinConstants = {
    O = {
      #Wpp
      -0.028
    }
    H = {
      # Wss
      -0.072
    }
  }
\end{verbatim}

\subsection{SpinOrbit}
\label{sec:dftbp.SpinOrbit}

If present, specifies that $L \cdot S$ coupling should be included for
a calculation. Currently spin unpolarised and non-collinear spin
polarisation are supported, but not collinear spin polarisation. For
every atomic species present in the calculation the spin-orbit
coupling constants, $\xi$, for that atom must be specified for all
shells present.  The constants must be ordered with respect to the
list of shells for the given atom.

In the case that the spin-orbit constant for an $s$ orbital has been
set to be a non-zero value the code prints a warning. For periodic
systems, use of this keyword automatically prevents the folding by
inversion normally used in \islcb{SupercellFolding}{sec:dftbp.SupercellFolding}, but manually
specified \is{KPointsAndWeights} should {\em not} be reduced by
inversion.

Example (GaAs):
\begin{verbatim}
  SpinOrbit = {
    Ga [eV] = {0.0 0.12 0.0} # s p d shells
    As [eV] = {0.0 0.32703} # s p shells
  }
\end{verbatim}

The additional option in this block, \is{Dual}, sets whether to use a block
population for the local spin matrices consistent with the dual populations of
Han~{\it et al.}~\cite{han-PRB-73-045110} or the conventional on-site part of
the single particle density matrix. The default value of this option is
\is{Yes}, also giving extra information regarding atomic orbital moments in the
detailed output.

\subsection{Solver}
\label{sec:dftbp.Solver}

Currently the following LAPACK~3.0~\cite{lapack3} eigensolver methods are always
available:
\begin{itemize}
\item \iscb{QR}\\ (QR decomposition based solver)
\item \iscb{DivideAndConquer}\\ (this requires about twice the memory of the
  other solvers)
\item \iscb{RelativelyRobust}\\ (using the subspace form but calculating all
  states)
\item \iscb{MAGMA}\\ (Only available for \dftbp{} binaries compiled with
  MAGMA~\cite{tdb10, tnld10, dghklty14} GPU support. \textbf{WARNING:} this is
  currently an experimental feature, so should be used with care.)
\end{itemize}
None of these solvers need any parameters or properties to be specified.

Example:\invparskip
\begin{verbatim}
  Solver = DivideAndConquer {}
\end{verbatim}

For ScaLAPACK enabled compilation, all three solvers are also available for MPI
parallel use.

If {\dftbp} is compiled with the ELSI library also included~\cite{YU2018267},
the additional \kw{ELPA}, \kw{OMM}, \kw{PEXSI} and \kw{NTPoly} solvers also
become available.

Note: The ELSI-solvers are not tested with multiple OpenMP-threads. Therefore,
{\dftbp} will stop with an error, if an ELSI-solver has been selected and the
maximal number of allowed threads is greater than one. (You can control the
number of allowed OpenMP-threads via the \verb|OMP_NUM_THREADS| environment
variable.)

\subsubsection{ELPA}

This is available with either single or two stage solution methods (the second
of these should be more efficiently parallel for large problems).

Example:\invparskip
\begin{verbatim}
  Solver = ELPA {
    Mode = 2
  }
\end{verbatim}

One caveat for this solver is that the number of parallel groups (see
p.~\pref{sec:dftbp.Parallel}) must match the number of k-points (times 2 in the
case of collinear spin polarisation). Calculations without k-points can use
either one or two groups in the case of collinear spin polarisation.

This solver can optionally use the ELSI internal sparse interface (\kw{Sparse} =
\is{Yes}), but this does not reduce memory usage and is mainly for testing
purposes.

\subsubsection{OMM}

This method minimises the single particle density matrix, so does not make band
structure information available. It is only stable for insulating grounds
states, i.e., systems with a HOMO-LUMO (band) gap.

The orbital minimisation method has four options:

\begin{ptable}
  \kw{nInterationsELPA} & i& & 5 & \\
  \kw{Tolerance} & r& & 1E-10 & \\
  \kw{Choleskii} & l& & Yes & \\
  \kw{Sparse}    & l& & No & \\
\end{ptable}
\begin{description}
\item[\is{nInterationsELPA}] Number of initial iterations to be performed with
  ELPA before the OMM method starts.
\item[\is{Tolerance}] Minimisation tolerance for this solver, larger values are
  faster by may be less stable.
\item[\is{Choleskii}] Whether the overlap is Choleskii factorised before
  applying OMM. This may increase stability of this method.
\item[\is{Sparse}] Whether the code should use the sparse matrix interface to
  ELSI solvers. This does not substantially improve memory usage in this case as
  internally the dense problem is solved with libOMM.
\end{description}

\subsubsection{PEXSI}

The PEXSI solver directly calculates the density matrix, so does not make band
structure information or Mermi free energy available. The scaling with system
size is better than the other solvers available in \dftbp{}, increasing as
$O(N_\mathrm{atom}^{d/2 + 1/2})$ where $d$ is the effective dimensionality of
the system. Hence for three dimensional structures it will scale as $O(N^2)$ for
general systems.
\begin{ptable}
  \kw{Method}       & i & (PEXSI $>$ 2.5) & 3 & \\
  \kw{Poles}        & i & \is{Method} = 3 & 30 & \\
  \kw{ProcsPerPole} & i & & 1  & \\
  \kw{muPoints}     & i & & 2  & \\
  \kw{SymbolicFactorProcs} & i & & 1  & \\
  \kw{SpectralRadius} & r & & 10 & \\
  \kw{Sparse} & l& & No & \\
  \kw{Threshold} & l& Sparse = No & 1E-15 & \\
\end{ptable}
\begin{description}
\item[\is{Method}] used for pole expansion. ELSI v2.5 only supports methods 1
  and 2, but v2.6 onwards offers method 3 as well (see the
  \href{https://wordpress.elsi-interchange.org/index.php/download/}{ELSI user
    manual} for details of the methods). Methods 1 and 3 can calculate the
  electron entropy, so enable the Mermin energy to be evaluated. When compiled
  with ELSI v2.5, the \dftbp{} default is \is{Method}=2, while for later
  versions method 3 is the default.
\item[\is{Poles}] number of poles for the complex plane calculation. \is{Method}=1
  defaults to 60 poles, 2 has 20 and method 3 uses 30 by default.
\item[\is{ProcsPerPole}] processors used to calculate the inversion at each pole.
\item[\is{muPoints}] number of processors used to search for the Fermi level.
\item[\is{SymbolicFactorProcs}] number of processors to use in evaluating the
  factorisation pattern of matrices.
\item[\is{SpectralRadius}]\modif{\modtype{Energy}} extension of the complex
  contour.
\item[\is{Sparse}] Whether the code should use the sparse ELSI matrix interface.
\item[\is{Threshold}] Sets the threshold to convert dense matrices to the
  internal sparse representation that ELSI uses. This may be useful in the case
  of matrix factorisation issues inside the solver.
\end{description}

\subsubsection{NTPoly}

This method constructs the single particle density matrix via a purification
method based on matrix polynomials (hence requires insulating systems). The
solver does not make band structure information available, but can be linear
scaling in both time and memory depending on settings and system. Currently the
solver does not support spin polarisation or k-points.

This solver has several options:
\begin{ptable}
  \kw{PurificationMethod} & i& & 2 & \\
  \kw{Tolerance} & r& & 1E-5 & \\
  \kw{Truncation} & r& & 1E-10 & \\
  \kw{Sparse} & l& & No & \\
  \kw{Threshold} & l& Sparse = No & 1E-15 & \\
\end{ptable}
\begin{description}
\item[\is{PurificationMethod}] Allowed choices are 0 for canonical purification,
  1 for trace correcting purification, 2 for 4$^\mathrm{th}$ order trace
  resetting purification, and 3 for generalised hole-particle canonical
  purification.
\item[\is{Tolerance}] Iterative convergence tolerance for this solver, larger
  values are faster by may be less stable.
\item[\is{Truncation}] Tolerance below which matrix elements in the density
  matrix are dropped to enforce sparsity.
\item[\is{Sparse}] Whether the code should use the sparse matrix ELSI interface.
\item[\is{Threshold}] Sets the threshold to convert dense matrices to the
  internal sparse representation that NTPoly uses.
\end{description}

The default choices of \is{Tolerance} and \is{Truncation} lead to an accurate,
but slow, solutions. Alternatively linear scaling can be achieved at smaller
system sizes with a larger choice of these values. Values in the range of 1E-3
and 1E-6 for \is{Tolerance} and \is{Truncation} may be suitable (but test the
quality of the solutions).

\subsection{Filling}
\label{sec:dftbp.Filling}

There are currently two types of filling supported (see below). Both have common
options:

\begin{ptable}
  \kw{Temperature} & r & AdaptFillingTemp = No & 0.0 & \\
  \kw{IndependentKFilling} & l& Periodic = Yes & No & \\
  \kw{FixedFermiLevel} & (1|2)r & & - & \\
\end{ptable}
\begin{description}
\item[\is{Temperature}]\modif{\modtype{energy}} Electron temperature in energy
  units. This property is ignored for thermostated MD runs, if the
  \is{AdaptFillingTemp} property of the thermostat has been set to \is{Yes} (See
  p.~\pref{sec:dftbp.Thermostat}).
\item[\is{IndependentKFilling}] Causes the occupation of the eigenstates to be independently
  determined for each $k$-point, thus preventing electron transfer between the $k$-points. Please
  note that the value for the Fermi level printed out by the code is meaningless in that case, since
  there is no common Fermi level for all $k$-points. This option is incompatible with use of the
  \is{FixedFermiLevel} keyword.
\item[\is{FixedFermiLevel}]\modif{\modtype{energy}} Can be used to fix the
  Fermi-level (total chemical potential, $\mu$) of the electrons in the
  system. For collinear spin polarisation, values for up and down spin channels
  are required. Otherwise only a single global chemical potential is
  required. If this option is present, the total charge and the total spin of
  the system are not conserved (settings in the options \is{Charge} and
  \is{UnpairedElectrons} will be ignored). If a fixed chemical potential is
  used, the output {\em force related energy} includes the contribution to the
  free energy, $- N \mu$, hence if differentiated will give the forces and
  stresses (if periodic).
\end{description}

\subsubsection{Fermi\cb}
\label{sec:dftbp.Fermi}

Fills the single particle levels according to a Fermi distribution. When using a
finite temperature, the Mermin free energy (which the code prints) should be
used instead of the total energy. This is given by $E - TS$, where the electron
entropy $S$ is used.

Example:
\invparskip
\begin{verbatim}
  Filling = Fermi {
    Temperature [K] = 300
  }
\end{verbatim}

\subsubsection{MethfesselPaxton\cb}
\label{sec:dftbp.MethfesselPaxton}

Produces a Fermi-like distribution but with much lower electron
entropy~\cite{methfessel-PRB-40-3616}. This is useful for systems that require
high electron temperatures (for example when calculating metallic systems). There
is an additional option for this type of filling:

\begin{ptable}
  \kw{Order} &i &  & 2 & \\
\end{ptable}
\begin{description}
\item[\is{Order}] Order of the Methessel-Paxton scheme, the order must be
  greater than zero, and the 1st order scheme is equivalent to Gaussian filling.
\end{description}

\textbf{Note:} Due to the non-monotonic behaviour of the Methfessel-Paxton filling function, the
position of the Fermi-level is not necessary unique for a given number of electrons. Therefore,
different fillings, band entropies, and Mermin free energies may result, depending which one has
been found by the Fermi-level search algorithm. The differences, however, are usually not physically
significant.

\subsection{Time-independent DFTB (TI-DFTB) excited states}

\subsubsection{$\boldsymbol\Delta$DFTB}
The time-independent approach $\Delta$DFTB modifies the SCC loop to allow
variational relaxation of the lowest singlet excited state of closed-shell
singlet species.\cite{irle-JCTC-12-313} The KS spin orbitals are self-constently
optimized under non-Aufbau orbital occupation constraints via promotion of an
electron from the highest occupied molecular orbital (HOMO) of the ground state
to lowest unoccupied molecular orbital (LUMO) at each SCC
iteration. Specifically, the \kwcb{NonAufbau} block in \dftbp manipulates the
electronic filling to emulate a singlet excited state. By modifying the
occupation pattern, the target excited state density ($\rho_e$) is can be
constructed according to the standard formula for obtaining the electron density
from the occupied colinear spin orbitals,
\begin{equation*}
  \rho_e(r) = \sum_\sigma \sum_{i \epsilon occ(\sigma)}\mid\phi_i^\sigma(r)\mid^2.
\end{equation*}
Because the constraint is placed on the spin orbitals rather than the MO's, this
method introduces heavy spin contamination to the system by the arbitrary
promotion of $\alpha$ spin electrons.  This issue is circumvented with spin
purification, where the triplet state energy is subtracted from twice the mixed
singlet state energy, in accordance with the Ziegler sum
rule,\cite{baerends-TCA-43-261}
\begin{equation*}
  E(S_1) = 2E[\{\phi_i^\sigma\}_m]-E[\{\phi_i^\sigma\}_t].
\end{equation*}
A non-spin-purified calculation can be performed if \is{SpinPurify = No}.
$\Delta$DFTB calculations can be more challenging to converge than their
ground-state analogs. To aid convergence, one can first run a ground-state
intial guess with \is{GroundGuess = Yes}.  The DIIS mixer is often an effective
alternative mixer for converging $\Delta$DFTB calculations in cases where the
calculation fails to converge with the Broyden mixer.

\label{sec:dftbp.NonAufbau}

This is a singlet excitation from a spin zero ground state system, so
\is{SpinPolarisation} should not be set in the input for this type of
calculation, but spin constants (p.~\pref{sec:dftbp.SpinConstants}) {\em are}
required. The calculation can be modified by including the following properties
in the \iscb{nonAufbau} block:
\begin{ptable}
  \kw{SpinPurify} & l & & Yes & \\
  \kw{GroundGuess} & l & & No & \\
\end{ptable}
\begin{description}
\item[\is{SpinPurify}] Calculates both a triplet and mixed state for
  spin-purification. If set to \is{No}, only the mixed state is calculated. The
  mixed state is formed by promoting an $\alpha$ electron from the HOMO to the
  LUMO.
\item[\is{GroundGuess}] Calculates the ground state energy prior to a non-Aufbau
  calculation.  The SCC loop will run three times if \is{NonAubau},
  \is{SpinPurify}, and \is{GroundGuess} are set to \is{Yes}, but only twice if
  \is{SpinPurify = No}.
\end{description}

The default settings, if \kwcb{NonAufbau} is included in the input, is
equivalent to:
\invparskip
\begin{verbatim}
  NonAufbau = {
    GroundGuess = No
    SpinPurify = Yes
  }
\end{verbatim}
This performs two self-consistent calculations at each geometry step, and then
evaluates the energy, forces and charges in the $\Delta$DFTB first excited
state. Note that this is a singlet state, hence only the {\em total} charge is
relevant in Mulliken populations.

If instead the ground state is also requested (\is{SpinPurify} = \is{Yes} by
default):
\invparskip
\begin{verbatim}
  NonAufbau = {
    GroundGuess = Yes
  }
\end{verbatim}
The \verb|detailed.out| will then contain energies of the ground, triplet and
the spin purified singlet excitation, along with estimates for the excitation
energies from the $S_o$ to the $T_1$ and $S_1$ states:
\begin{verbatim}
S0 -> T1:                            0.2104101411 H            5.7256 eV
S0 -> S1:                            0.2615036445 H            7.1159 eV
\end{verbatim}

\subsection{SlaterKosterFiles}
\label{sec:dftbp.SlaterKosterFiles}

There are two different ways to specify the Slater-Koster files for
the atom type pairs, explicit specification and using the
\iscb{Type2FileNames} method.

\subsubsection{Explicit specification}

Every pairwise permutation atomic types, connected by a dash, must
occur as a property with the name of the corresponding file as an
assigned value.

Example (GaAs):
  \invparskip
\begin{verbatim}
  SlaterKosterFiles = {
    Ga-Ga = "./Ga-Ga.skf"
    Ga-As = "./Ga-As.skf"
    As-Ga = "./As-Ga.skf"
    As-As = "./As-As.skf"
  }
\end{verbatim}

If you treat shells from different species as shells of one atom by
using the \iscb{SelectedShells} keyword in the
\iscb{MaxAngularMomentum} block, you have to specify more than one
file name for certain species pairs. (For details see the description
about the \iscb{MaxAngularMomentum} keyword.)

\subsubsection{Type2FileNames\cb}
\label{sec:dftbp.Type2FileNames}

You can use this method to generate the name of the Slater-Koster
files automatically using the element names from the input
geometry. You have to specify the following properties:
\begin{ptable}
  \kw{Prefix} & s &  & "" & \\
  \kw{Separator} &s &  & "" & \\
  \kw{Suffix} & s & & "" & \\
  \kw{LowerCaseTypeName} & l & & No & \\
\end{ptable}
\begin{description}
\item[\is{Prefix}] Prefix before the first type name, usually the path.
\item[\is{Separator}] Separator between the type names.
\item[\is{Suffix}] Suffix after the name of the second type, usually
  extension.
\item[\is{LowerCaseTypeName}] If the name of the types should be
  converted to lower case. Otherwise they are used in the same way, as
  they were specified in the geometry input.
\end{description}

Example (for producing the same file names as in the previous section):
\invparskip
\begin{verbatim}
  SlaterKosterFiles = Type2FileNames {
    Prefix = "./"
    Separator = "-"
    Suffix = ".skf"
    LowerCaseTypeName = No
  }
\end{verbatim}

The \is{Type2FileNames} method can not be used if an extended basis
was defined with the \is{SelectedShells} method.


\subsection{KPointsAndWeights}
\label{sec:dftbp.KPointsAndWeights}
\index{Brillouin-zone sampling}

The $k$-points for the Brillouin-zone integration can either be specified
explicitly, or automatically for supercells by using the \iscb{SupercellFolding}
or \iscb{KLines} methods for supercells or either \iscb{HelicalUniform} or
\iscb{HelicalSampled} for helical boundary conditions.

\textbf{ Note:} For the automatic grid methods, the \is{KPointsAndWeights}
keyword is not allowed to have a modifier.

\subsubsection{Explicit specification}

Each $k$-point and its weight in the integral should be specified, for
supercells this requires that four real numbers must be specified for
each point: The coordinates of the given $k$-point followed by its
weight, while for helical coordinates there are two coordinates (along
the helical axis and with respect to the rotation around the axis) and
the weight of the point. By default, coordinates are specified in
fractions of the reciprocal lattice vectors. If the modifier
\is{absolute} is set for the \is{KPointsAndWeights} keyword, absolute
$k$-point coordinates in atomic units are instead expected.  The sum
of the k-point weights is automatically normalised by the program.
\begin{verbatim}
  KPointsAndWeights = {   # 2x2x2 MP-scheme
    0.25  0.25  0.25    1.0
    0.25  0.25 -0.25    1.0
    0.25 -0.25  0.25    1.0
    0.25 -0.25 -0.25    1.0
  }
\end{verbatim}

\subsubsection{SupercellFolding\cb}
\label{sec:dftbp.SupercellFolding}

This method generates a sampling set containing all the special
k-points in the Brillouin zone related to points that would occur in
an enlarged supercell repeating of the current unit cell.  If two
$k$-points in the BZ are related by inversion, only one (with double
weight) is used (in the absence of spin-orbit coupling this is
permitted by time reversal symmetry). The \iscb{SupercellFolding}
method expects 9 integers and 3 real values as parameters:
\begin{equation*}
  \begin{array}{ccc}
    n_{11} & n_{12} & n_{13} \\
    n_{21} & n_{22} & n_{23} \\
    n_{31} & n_{32} & n_{33} \\
    s_{1} & s_{2}   & s_{3} \\
  \end{array}
\end{equation*}
The integers $n_{ij}$ specify the coefficients used to build the
supercell vectors $\mathbf{A}_i$ from the original lattice vectors
$\mathbf{a}_j$:
\begin{equation*}
  \mathbf{A}_i = \sum_{j=1}^3 n_{ij}\, \mathbf{a}_j
  \text.
\end{equation*}
The real values, $s_i$, specify the point in the Brillouin-zone of the
super lattice, in which the folding should occur. The coordinates must
be given in relative coordinates, in the units of the reciprocal
lattice vectors of the super lattice.

The original $l_1\times l_2\times l_3$ Monkhorst-Pack
sampling\index{Monkhorst-Pack scheme} \cite{monkhorst-prb-13-5188} for
cubic lattices corresponds to a uniform extension of the lattice:
\begin{equation*}
  \begin{array}{ccc}
    l_1 & 0 & 0\\
    0 & l_2 & 0 \\
    0 & 0& l_3 \\
    s_1 & s_2 & s_3
  \end{array}
\end{equation*}
where $s_i$ is $0.0$, if $l_i$ is odd, and $s_i$ is $0.5$ if $l_i$ is
even. For the $2\times2\times3$ scheme, you would write for example
\begin{verbatim}
  # 2x2x3 MP-scheme according original paper
  KPointsAndWeights = SupercellFolding {
     2    0    0
     0    2    0
     0    0    3
     0.5  0.5  0.0
  }
\end{verbatim}

To use k-points for hexagonal lattices which are consistent with the
erratum to the original paper \cite{monkhorst-prb-16-1748}, you should
set the shift for the unique ``$c$'' direction, $s_3$, in the same way
as in the original scheme. The $s_1$ and $s_2$ shifts should be set to
be $0.0$ independent of whether $l_1$ and $l_2$ are even or odd.  So,
for a $2\times3\times4$ sampling you would have to set
\begin{verbatim}
  # 2x3x4 MP-scheme according modified MP scheme
  KPointsAndWeights = SupercellFolding {
     2    0    0
     0    3    0
     0    0    4
     0.0  0.0  0.5
  }
\end{verbatim}

It is important to note that \dftbp{} does not take the symmetry of
your system explicitly into account. For small high symmetric systems
with a low number of $k$-points in the sampling this could eventually
lead to unphysical results. (Components of tensor
properties--e.g.\ forces--could be finite, even if they must vanish
due to symmetry reasons.) For those cases, you should explicitly
specify $k$-points with the correct symmetry.


\subsubsection{KLines\cb}
\label{sec:dftbp.KLines}

This method specifies $k$-points lying along arbitrary lines in the
Brillouin zone. This is useful when calculating the \index{band
  structure calculation}{band structure} for a periodic system. (In
that case, the charges should be initialised from the saved charges of
a previous calculation with a proper $k$-sampling. Additionally for
SCC calculations the number of SCC cycles should be set to 1, so that
only one diagonalisation is done using the initial charges.)

The \iscb{KLines} method accepts for each line an integer specifying
the number of points along the line segment, and 3 real values
specifying the end point of the line segment. The line segments do not
include their starting points but their end points. The starting point
for the first line segment can be set by specifying a (zeroth) segment
with only one point and with the desired starting point as end point.
The unit of the $k$-points is determined by any modifier of the
\is{KPointsAndWeights} property. (Default is relative coordinates.)

Example:
\invparskip
\begin{verbatim}
  KPointsAndWeights [relative] = KLines {
    1   0.5  0.0  0.0    # Setting (and calculating) starting point 0.5 0.0 0.0
   10   0.0  0.0  0.0    # 10 points from 0.5 0.0 0.0  to  0.0 0.0 0.0
   10   0.5  0.5  0.5    # 10 points from 0.0 0.0 0.0 to 0.5 0.5 0.5
    1   0.0  0.0  0.0    # Setting (and calculating) a new starting point
   10   0.5  0.5  0.0    # 10 points from 0.0 0.0 0.0 to 0.5 0.5 0.0
  }
\end{verbatim}

\textbf{Note:} Since this set of k-points probably does not correctly integrate the Brillouin zone,
the default value of \kw{MaxSccIterations} is set to be 1 in this case.

\htcbsubsubsection{HelicalUniform}

This method specifies $k$-points lying along the generalized reciprocal vector
of the Brillouin zone of a helical cell and around the order-$n$ rotational axis
(currently the $k$-points that exactly represent the $C_n$ rotation are
used). The \iscb{HelicalUniform} method expects 1 integer and 1 real value as
parameters, where the first value specifies the number of sampling points along
the helical axis, while the second gives the shift (analogous to the three
dimensional case of \iscb{SupercellFolding}) in this direction. A shift of $0.5$
appears to give more rapid convergence of the grid.

Example:
\invparskip
\begin{verbatim}
  KPointsAndWeights = HelicalUniform {80 0.5}
\end{verbatim}

\htcbsubsubsection{HelicalSampled}

Instead of exactly integrating around the $C_n$ rotation in $k$-space, the
\iscb{HelicalSampled} method allows for a sampled integration. It expects 4
values: the first two are the number of sample points along the helix and around
the rotation respectively, while the second two are shifts in the grid.

Example:
\invparskip
\begin{verbatim}
  KPointsAndWeights = HelicalSampled {20 4 0.5 0.25}
\end{verbatim}

There are several things to note here: firstly, the second grid (values 4 and
0.25 in the example above) approximates the integration around the $C_n$
symmetry, so should only be used where the order of this rotation axis is
large. Secondly, non-zero shifts in the grid, particularly for small numbers of
sampling points, are likely to be unphysical as are shifts for the grids at the
exact order of the rotation operation.

\subsection{OrbitalPotential}
\label{sec:dftbp.OrbitalPotential}

\index{DFTB+U}
\label{sec:DFTB+U}

Currently the \is{FLL} (fully localised limit) and
\is{pSIC}~\cite{hourahine07} (pseudo self interaction correction )
forms of the LDA+U corrections~\cite{petukhov-PRB-67-153106} are
implemented. These potentials effect the energy of states on
designated shells of particular atoms, usually increasing the
localisation of states at these sites. The \is{FLL} potential lowers
the energy of occupied states localised on the specified atomic shells
while raising the energy of unoccupied states. The the \is{pSIC}
potential corrects the local part of the self-interaction error and so
lowers the energy of occupied states (see Ref.~\cite{hourahine07} for
a discussion of the relation between these two potentials, and
possible choices for the UJ constant).  These particular corrections
are most useful for lanthanide/actinide $f$ states and some localised
$d$ states of transition metals (Ni$3d$ for example).

The \is{Functional} option chooses which correction to apply, followed
by a list of the specific corrections, listed as an atomic species and
the shells on that atom followed by the $U-J$ constant for that block
of shells.

\begin{verbatim}
  OrbitalPotential = {
   Functional = {FLL}
   Si = {
     Shells = {1 2} # sp block on the atom
     UJ = 0.124
   }
  }
\end{verbatim}

\subsection{ElectricField}
\label{sec:dftbp.ElectricField}

This tag contains the specification for an external electric
field. Electric fields can only be specified for SCC calculations. You
can apply the electric field of point charges\footnote{Only in
  calculations with fixed lattice constants.} and/or a homogeneous
external field (which may change harmonically in time). The
\is{ElectricField} block can currently contain either one or more
\kw{PointCharges} blocks and potentially an \kw{External} block.

\subsubsection{PointCharges}
\label{sec:dftbp.PointCharges}
The specification for \kw{PointCharges} has the following properties:
\begin{ptable}
  \kw{CoordsAndCharges} & (4r)+ & & - & \\
  \kw{GaussianBlurWidth} & r & Periodic = No & 0.0 & \\
\end{ptable}
\begin{description}
\item[\is{CoordsAndCharges}]\modif{\modtype{length}} Contains the
  coordinates and the charge for each point charge (four real values
  per point charge). A length modifier can be used to alter the units
  of the coordinates. The charge must be specified in proton
  charges. (The charge of an electron is -1.)

  If you read in a huge number of external charges the parsing time to
  process this data could be unreasonably long. You can avoid this by
  including the coordinates and the charges directly from an external
  file via the \kwcb{DirectRead} method (see the example in the next
  paragraph). Please note that when using this method the program will
  only read the specified number of records from the external file,
  and ignores any additional data (so do not leave comments in the
  external file for example). The external file should contain only
  one record (3 coordinates and 1 charge) per line.

\item[\is{GaussianBlurWidth}]\modif{\modtype{length}} Specifies the
  half width $\sigma$ of the Gaussian charge distribution, which is
  used to delocalise the point charges.  The energy of the coulombic
  interaction $E_{\text{C}}$ between the delocalised point charge $M$
  with charge $Q_M$ and the atom $A$ with charge $q_A$ is weighted by
  the error function as
  \begin{equation*}
    E_{\text{C}}(A,M) = \frac{q_A Q_M}{r_{AM}}\,
    \operatorname{erf}\left[\frac{r_{AM}}{\sigma}\right]
    \text,
  \end{equation*}
  where $r_{AM}$ is the distance between the point charge and the
  atom.

  A length modifier can be used to specify the unit for $\sigma$.

  Example:\invparskip
\begin{verbatim}
  ElectricField = {
    # 1st group of charges, with Gaussian delocalisation
    # We have 100000 charges, therefore we choose the fast reading method.
    PointCharges = {
      GaussianBlurWidth [Angstrom] = 3.0
      CoordsAndCharges [Angstrom] = DirectRead {
        Records = 100000
        File = "charges.dat"
      }
    }
    # 2nd group of charges, no delocalisation (sigma = 0.0)
    PointCharges = {
      CoordsAndCharges [Angstrom] = {
        3.3  -1.2  0.9      9.2
        1.2  -3.4  5.6     -3.3
      }
    }
  }
\end{verbatim}
\end{description}

\subsubsection{External}
\label{sec:dftbp.External}

Specifies a homogeneous external electric field. In the case of {\em
  periodic} calculations, a saw-tooth potential is currently used,
hence it is up to the user to guarantee that there is a vacuum region
isolating periodic copies of the system along the applied field
direction.  We suggest that you place the structure in the `middle' of
the unit cell if possible, to reduce the chances of atoms approaching
cell boundaries along the direction of the applied electric field. The
code will halt if atoms interact with periodic images of the unit cell
along the direction of the electric field.

The \kw{External} field keyword has the following options
\begin{ptable}
  \kw{Strength}  & r  &                                  & -   & \\
  \kw{Direction} & 3r &                                  &     & \\
  \kw{Frequency} & r  & \textrm{molecular dynamics used} & 0.0 & \\
  \kw{Phase}     & i  & Geometry step offset             & 0   & \\
\end{ptable}
\begin{description}
\item[\is{Strength}]\modif{\modtype{Electric field strength}}
  Specified strength of the applied field.
\item[\is{Direction}] Vector direction of the applied field (the code
  normalises this vector). In the case of periodic calculations,
  currently the system {\em must not} be continuous in this
  direction (see above).
\item[\is{Frequency}]\modif{\modtype{Frequency}} If using molecular
  dynamics, the field can be time varying with this frequency.
\item[\is{Phase}] Initial field phase in units of geometry steps, this is needed
  if restarting an MD run in an external field to give the offset in phase of
  the field after the specified number of steps from the old calculation. The
  applied field is of the form
  \begin{equation*}
    \mathbf{E}_0 \sin( \omega \Delta t (step + phase) )
  \end{equation*}
  where $\mathbf{E}_0$ is the field vector specified by \kw{Strength} and
  \kw{Direction}, $\omega$ is the angular \kw{Frequency} and $step$ is the
  current MD-step in the simulation, using the MD \kw{TimeStep} of $\Delta t$
  (see section \ref{sec:dftbp.VelocityVerlet}).
\end{description}

\subsection{Dispersion}
\label{sec:dftbp.Dispersion}

The \is{Dispersion} block controls whether DFTB interactions should be
empirically corrected for van der Waals interactions, since DFTB (and
SCC-DFTB) does not include these effects. Currently, six different
dispersion correction schemes are implemented (for the detailed
description of the methods see the following subsections):
\begin{itemize}
\item \is{LennardJones}: Dispersion is included via a Lennard-Jones
  potential between each pair of atoms. The parameters for the
  potential can either be entered by the user or the program can
  automatically take the parameters from the Universal Force Field
  (UFF)~\cite{rappe-JACS-114-10024}.
\item \is{SlaterKirkwood}: The dispersion interaction between atoms is
  taken from a Slater-Kirkwood polarisable atomic
  model~\cite{elstner-jcp-114-5149}.
\item \is{DftD3}: Dispersion is calculated as in the dftd3
  code~\cite{grimme-jcp-132-154104,grimme-jcp-32-1456-1465} (see
  section~\ref{sec:dftbp.DftD3}). Modification hydrogen bond interaction
  strengths (see section \ref{sec:dftbp.hcorr}).
\item \is{DftD4}: Dispersion is calculated using the D4
  model~\cite{caldweyher-jcp-147-034112,caldeweyher-jcp-150-154122} (see
  section~\ref{sec:dftbp.DftD4}).
\item \is{Ts}: Dispersion is calculated using the Tkatchenko-Scheffler model~\cite{TkatchenkoPRL09}
adapted for DFTB~\cite{StohrJCP16} (see section~\ref{sec:dftbp.DispTS}).
\item \is{Mbd} Dispersion is calculated using the MBD@rsSCS model~\cite{AmbrosettiJCP14}
adapted for DFTB~\cite{StohrJCP16} (see section~\ref{sec:dftbp.DispMBD}).
\end{itemize}

\subsubsection{LennardJones}
\label{sec:dftbp.LennardJones}

The Lennard-Jones dispersion model in \dftbp{} follows the method of
Ref.~\cite{zhechkov-JCTC-1-841}, using the following potential:
\begin{eqnarray*}
U_{ij}(r)&=&d_{ij}\left[-2\left(\frac{r_{ij}}{r}\right)^6 +
  \left(\frac{r_{ij}}{r}\right)^{12}\right]\qquad r >= r_0\\
U_{ij}(r)&=&U_0 + U_1 r^5 + U_2 r^{10}\qquad r < r_0\\
\end{eqnarray*}
where $r_0$ is the distance at which the potential turns from
repulsive to attractive. The parameters $d_{ij}$ and $r_{ij}$ are
built from atomic parameters $d_i$, $d_j$ and $r_i$, $r_j$ via the
geometrical mean ($d_{ij} = \sqrt{d_id_j}$,
$r_{ij}=\sqrt{r_ir_j}$). The parameters $U_0$, $U_1$, $U_2$ ensure a
smooth functional form at $r_0$.

The parameters $r_i$ and $d_i$ can either be taken from the parameters
of the UFF~\cite{rappe-JACS-114-10024} (as in
Ref.~\cite{zhechkov-JCTC-1-841}) or can be specified manually for each
species.

Example using UFF parameters:\invparskip
\begin{verbatim}
  Dispersion = LennardJones {
    Parameters = UFFParameters {}
  }
\end{verbatim}

Example using manually specified parameters:\invparskip
\begin{verbatim}
  Dispersion = LennardJones {
    Parameters {
     H {
       Distance [AA] = 2.886
       Energy [kcal/mol] = 0.044
     }
     O {
       Distance [AA] = 3.500
       Energy [kcal/mol] = 0.060
     }
    }
  }
\end{verbatim}

The UFF provides dispersion parameters for nearly every element of the
periodic table, therefore it can be used for almost all systems
``out of the box''. The parameters are also independent of the atomic
coordination number, allowing straight forward geometry relaxation or
molecular dynamics (unlike the current implementation of
Slater-Kirkwood type dispersion).


\subsubsection{SlaterKirkwood}
\label{sec:dftbp.SlaterKirkwood}
\label{sec:SlaterKirkwood}

A Slater-Kirkwood type dispersion model is also implemented in
\dftbp{} as described in
Ref.~\cite{elstner-jcp-114-5149}.\footnote{Please note, that
  Ref.~\cite{elstner-jcp-114-5149} contains two typos: equation (7)
  should read $C_6^{\alpha\beta} = \frac{2 C_6^\alpha C_6^\beta
    p_\alpha p_\beta}{p_\alpha^2 C_6^\beta + p_\beta^2 C_6^\alpha}$,
  in equation (9) the contribution from the dispersion should be
  $E_{\text{dis}} = -\frac{1}{2} \sum_{\alpha\beta}
  f(R_{\alpha\beta})C_6^{\alpha\beta}(R_{\alpha\beta})^{-6}$.}
This model requires atomic polarisation
values, van der Waals radii and effective charges for every atom in
your system. These parameters are dependent on the coordination of
each atom, hence values for different atoms of the same species may
vary depending on local environment.  You can supply these parameters
for the atoms in either of two ways, both using the
\is{PolarRadiusCharge} tag.

The first option is to specify the values within the
\is{PolarRadiusCharge} environment by providing three real values
(polarisability, van der Waals radius, effective charge) for each atom
separately.

Example:\invparskip
\begin{verbatim}
  Dispersion = SlaterKirkwood {
    # Using Angstrom^3 for volume, Angstrom for length and default
    # unit for charge (note the two separating commas between the units)
    PolarRadiusCharge [Angstrom^3,Angstrom,] = {
      # Polar      Radius     Chrg
      0.560000    3.800000    3.150000      # Atom 1: O
      0.386000    3.500000    0.800000      # Atom 2: H
      0.386000    3.500000    0.800000      # Atom 3: H
    }
  }
\end{verbatim}

Alternatively you can provide values for each atomic species in your
system, but must supply different values for different coordination
numbers using the \kwcb{HybridDependentPol} keyword. The code needs
specific parameters for each type of atom in environments with 0, 1,
2, 3, 4 or $\geqslant$5 neighbours. \dftbp{} then picks the
appropriate values for each atom based on their coordination in the
\emph{starting} geometry.  Two atoms are considered to be neighbours
if their distance is less than the sum of their covalent radii, hence
you need to supply the covalent radii for each atomic species using
the \kw{CovalentRadius} keyword. This is then followed by a
\kw{HybridPolarisations} block containing a list of six values for
atomic polarisabilities then six van der Waals radii and finally a
single hybridisation independent effective charge for that atomic
species.

Example:\invparskip
\begin{verbatim}
  Dispersion = SlaterKirkwood {
    PolarRadiusCharge = HybridDependentPol {
      O = {
        CovalentRadius [Angstrom] = 0.8
        HybridPolarisations [Angstrom^3,Angstrom,] = {
          # Atomic polarisabilities 0-5        van der Waals radii 0-5  chrg
          0.560 0.560 0.560 0.560 0.560 0.560  3.8 3.8 3.8 3.8 3.8 3.8  3.15
        }
      }
      H = {
        CovalentRadius [Angstrom] = 0.4
        HybridPolarisations [Angstrom^3,Angstrom,] = {
          # Atomic polarisabilities 0-5        van der Waals radii 0-5  chrg
          0.386 0.396 0.400 0.410 0.410 0.410   3.5 3.5 3.5 3.5 3.5 3.5   0.8
        }
      }
    }
  }
\end{verbatim}

\textbf{Warning:} For both methods of specifying the Slater-Kirkwood dispersion model the code keeps
the dispersion parameters fixed for each atom during the entire calculation. Even if the geometry
(and therefore the hybridisation) of atoms changes significantly during the calculation, the
parameters are unchanged. Therefore if atoms are able to move during your calculation (geometry
relaxation or molecular dynamics) you should \emph{always} check whether the coordination of your
atoms has changed during the run.

Furthermore, when using the \iscb{HybridDependentPol} method we
suggest that you first set the \is{StopAfterParsing} keyword in the
\is{ParserOptions} block to \is{Yes} (see p.~\pref{sec:dftbp.ParserOptions}) and
inspect the generated polarisabilities, radii and charges for every
atom in the \verb|dftb_pin.hsd| file. If fine tuning of the generated
values turns out to be necessary, you should replace the hybrid
dependent specification in the input file with corrected atom
specific values based on \verb|dftb_pin.hsd|.

In order to find suitable parameters for the Slater-Kirkwood model,
you should consult Ref.~\cite{elstner-jcp-114-5149} and further
references therein. Appendix \ref{app:dispconsts} contains values
which have already been used by some DFTB-users for a few elements.

\subsubsection{DftD3}
\label{sec:dftbp.DftD3}

The DFT-D3 dispersion correction in \dftbp{} is an implementation of the method
used in the code 'dftd3' by Stefan Grimme and coworkers.  It is based on the
{\it ansatz} described in Refs.~\cite{grimme-jcp-132-154104} and
\cite{grimme-jcp-32-1456-1465}.

\textbf{Note:} the \dftbp{} binary must be compiled with the DFT-D3 library
enabled to use this feature.

This dispersion correction for DFTB adds a contribution to the general
Kohn-Sham-like energy
\begin{equation*}
  E_{\text{DFTB-D3}} = E_{\text{DFTB}} + E_{\text{disp}}
\end{equation*}
with $E_{\text{DFTB}}$ being the DFTB total energy and $E_{\text{disp}}$ the
dispersion energy. The latter contains two-body and optional three-body
contributions:
\begin{equation*}
  E_{\text{disp}} = E_{\text{disp}}^{(2)} + E_{\text{disp}}^{(3)}.
\end{equation*}

The form of the two-body contribution can change depending on the chosen damping
factor:
\begin{itemize}
\item Becke-Johnson damping function:
  \begin{equation*}
    E_{\text{disp}}^{(2)} = -\frac{1}{2} \sum_{A\neq B} \sum_{n=6,8} s_n
    \frac{C_n^{AB}}{r_{AB}^n + f(R_0^{AB})}
  \end{equation*}
  with
  \begin{equation*}
    f(R_0^{AB}) = a_1 R_0^{AB} + a_2 \text.
  \end{equation*}

\item Zero-damping (dispersion at short distances is damped to zero):
  \begin{equation*}
    E_{\text{disp}}^{(2)} = -\frac{1}{2} \sum_{A \neq B} s_n
    \frac{C_n^{AB}}{r_{AB}^n} f_{d,n}(r_{AB})
  \end{equation*}
  with
  \begin{equation*}
    f_{d,n} = \frac{1}{1 + 6(r_{AB}/(s_{\text{r},n} R_0^{AB}))^{-\alpha_n}}
  \end{equation*}
\end{itemize}

In order to adjust the dispersion for various energy functionals, the choice of
$s_6$, $s_8$ and the damping parameters $a_1$ and $a_2$ (for
Becke-Johnson-damping) or $s_{\text{r},6}$ and $\alpha_6$ (for zero damping) are
treated as functional-dependent values. All other parameters are fixed based on
these parameters.

As the DFTB energy functional is largely determined by the underlying
parameterisation (the Slater-Koster-files) and the chosen DFTB model (e.g.\
non-scc, scc, 3rd order, etc.), there are no universal parameter choices which
can be used with all settings, but some relevant choices for various
parameterisation are given in Appendix \ref{app:dftd3const}.

\textbf{Note:} for the version 6 or earlier of the \dftbp{} input parser (see
section~\ref{sec:dftbp.ParserOptions}) the default values of these parameters
are set to be appropriate for DFTB3. But from parser version 7 onwards, no
default values are set.

Example using adjusted parameters with Becke-Johnson damping:
\begin{verbatim}
  Dispersion = DftD3 {
    Damping = BeckeJohnson {
      a1 = 0.5719
      a2 = 3.6017
    }
    s6 = 1.0
    s8 = 0.5883
  }
\end{verbatim}

Example using zero-damping:
\begin{verbatim}
  Dispersion = DftD3 {
    Damping = ZeroDamping {
      sr6 = 0.7461
      alpha6 = 14.0
    }
    s6 = 1.0
    s8 = 3.209
  }
\end{verbatim}

\subsubsection{DftD3 optional settings}
Apart from the functional dependent dispersion parameters, you can also adjust
the additional parameters as shown below. The default values for these
parameters are taken to be the same as in the dftd3 code.

\begin{ptable}
  \kw{Cutoff} & r & & $\sqrt{9000}$ & \\
  \kw{CutoffCN} & r & & $40$ & \\
  \kw{Threebody} & l & & No & \\
  \kw{HHRepulsion} & l & & No & \\
  \kw{AtomicNumbers} & m & & \cb & \\
\end{ptable}
\begin{description}
\item[\is{Cutoff}] \modif{\modtype{length}} Cutoff distance when calculating
  two-body interactions.

\item[\is{CutoffCN}] \modif{\modtype{length}} Cutoff distance when calculating
  three-body interactions.

\item[\is{Threebody}] Whether three-body contributions should be included in the
  dispersion interactions.

\item[\is{AtomicNumbers}]
  Allows overwrite of atomic numbers associated with a species name.
  \begin{verbatim}
  # Default species (H-Og) are implied
  AtomicNumbers {
    D = 1  # set deuterium as hydrogen
    S = 6  # overwrite sulfur as carbon
  }
  \end{verbatim}

\item[\is{HHRepulsion}] Required when calculating the
  DFTB3-D3H5~\cite{rezac-jctc-13-2017} modification to D3 dispersion (see
  section~\ref{sec:dftbp.hcorr} for details and parameter values). This keyword
  enables an additional short range repulsion term in all hydrogen--hydrogen
  pairs~\cite{rezac-jctc-8-2012} which prevents them from approaching too
  closely together.
\end{description}

\subsubsection{DftD4}
\label{sec:dftbp.DftD4}

The DFT-D4 dispersion correction in \dftbp{} is an implementation of the D4
model by Stefan Grimme and coworkers. It is based on the method described in
Ref.~\cite{caldeweyher-jcp-150-154122}.

This dispersion correction for DFTB adds a contribution to the general
Kohn--Sham-like energy,
\begin{equation*}
  E_{\text{DFTB-D4}} = E_{\text{DFTB}} + E_{\text{disp}},
\end{equation*}
with $E_{\text{DFTB}}$ being the DFTB total energy and $E_{\text{disp}}$ the
dispersion energy. The latter contains two-body and optional three-body
contributions:
\begin{equation*}
  E_{\text{disp}} = E_{\text{disp}}^{(2)} + E_{\text{disp}}^{(3)}.
\end{equation*}

The D4 model uses the Becke--Johnson damping function for two-body contributions:
\begin{equation*}
   E_{\text{disp}}^{(2)} = -\frac{1}{2} \sum_{A\neq B} \sum_{n=6,8,10} s_n
   \frac{C_n^{AB}}{r_{AB}^n + f(R_0^{AB})}
\end{equation*}
with
\begin{equation*}
   f(R_0^{AB}) = a_1 R_0^{AB} + a_2.
\end{equation*}

The zero-damping function for three-body contributions is:
\begin{equation*}
   E_{\text{disp}}^{(3)} = -\sum_{A}\sum_{\substack{B \\B < A}}\sum_{\substack{C\\ C < B}} s_9
   \frac{  \left(3\cos\theta_A\cos\theta_B\cos\theta_C+1\right)\sqrt{C_6^{AB}C_6^{BC}C_6^{CA}}}{(r_{AB}r_{BC}r_{CA})^3} f^{(3)}(r_{AB})
\end{equation*}
with
\begin{equation*}
   f^{(3)} = \frac{1}{1 + 6\left(\frac{f(R_0^{AB})f(R_0^{BC})f(R_0^{CA})}{r_{AB}r_{BC}r_{CA}}\right)^{\alpha/3}}.
\end{equation*}

In order to adjust the dispersion for various energy functionals, the choice of
$s_8$ and the damping parameters $a_1$ and $a_2$ are
treated as functional-dependent values. All other parameters are fixed based on
these parameters.
Depending on the choice of the $s_9$ parameter non-additive triple-dipole
contributions will be evaluated. Including non-additive effects usually improves
the description of dispersion interactions, but is also more expensive.

As the DFTB energy functional is largely determined by the underlying
parameterisation (the Slater--Koster-files) and the chosen DFTB model (e.g.\
non-scc, scc, 3rd order, etc.), there are no universal parameter choices which
can be used with all settings, but some relevant choices for various
parameterisation are given in Appendix \ref{app:dftd4const}.

\subsubsection{DftD4 settings}
Beside the functional dependent dispersion parameters, the options shown below
can be adjusted in the input.

\begin{ptable}
  \kw{s6} & r & & 1.0 & \\
  \kw{s8} & r & & & \\
  \kw{s10} & r & & 0.0 & \\
  \kw{s9} & r & & & \\
  \kw{a1} & r & & & \\
  \kw{a2} & r & & & \\
  \kw{alpha} & r & & 16.0 & \\
  \kw{WeightingFactor} & r & & 6.0 & \\
  \kw{ChargeSteepness} & r & & 2.0 & \\
  \kw{ChargeScale} & r & & 3.0 & \\
  \kw{CutoffInter} & r & & 64 & \\
  \kw{CutoffThree} & r & & 40 & \\
  \kw{CoordinationNumber} & m & & \is{Cov} & \pref{sec:dftbp.CoordinationNumber} \\
  \kw{ChargeModel} & m & & \is{EEQ} & \pref{sec:dftbp.ChargeModel} \\
  \kw{AtomicNumbers} & m & & \cb & \\
\end{ptable}

\begin{description}
\item[\is{s8}, \is{s9}, \is{a1}, \is{a2}] Functional dependent dispersion
  parameters, see equations above.

\item[\is{s6}] Parameter for scaling pairwise dipole--dipole dispersion
  interaction, should always be set to 1.0.

\item[\is{s10}] Parameter for pairwise quadrupole--quadrupole dispersion
  interactions, should be kept set to 0.0.

\item[\is{alpha}] Zero damping exponent for three-body damping function, default
  is 16 as in DFT-D3.

\item[\is{WeightingFactor}] Coordination number based interpolation, 4.0 used in
  DFT-D3, default for D4 is 6.0.

\item[\is{ChargeScale}] Maximum possible charge scaling, used as exponential
  value, should be kept set to 3.0.

\item[\is{ChargeSteepness}] Steepness of the charge scaling function, should be
  kept set to 2.0.

\item[\is{CutoffInter}] \modif{\modtype{length}} Cutoff distance when
  calculating two-body interactions.

\item[\is{CutoffThree}] \modif{\modtype{length}} Cutoff distance when
  calculating three-body interactions.

\item[\is{AtomicNumbers}]
  Allows overwrite of atomic numbers associated with a species name.

\end{description}

\subsubsection{DftD4 ChargeModel}
\label{sec:dftbp.ChargeModel}

This implementation of DFT-D4 supports only the \is{EEQ\cb} method to initialize
the charge model with an electronegativity equilibration~(EEQ) model.\cite{rappe1991}
For each species four parameters (\kw{Chi}, \kw{Gam}, \kw{Kcn}, and \kw{Rad})
have to be supplied in a \is{Values\cb} method, since the model
is instanciated inside the \is{DftD4\cb} method, \is{Defaults\cb} for all elements
up to 86 can be supplied automatically.\cite{caldeweyher-jcp-150-154122}

\begin{ptable}
  \kw{Chi} & m & & \is{Defaults} & \\
  \kw{Gam} & m & & \is{Defaults} & \\
  \kw{Kcn} & m & & \is{Defaults} & \\
  \kw{Rad} & m & & \is{Defaults} & \\
  \kw{Cutoff} & r & & 40 & \\
  \kw{EwaldParameter} & r & & 0.0 & \\
  \kw{EwaldTolerance} & r & & 1.0e-9 & \\
  \kw{CoordinationNumber} & m & & \is{Erf} & \pref{sec:dftbp.CoordinationNumber} \\
\end{ptable}

\begin{description}

\item[\is{Chi}] Electronegativities of all species.

\item[\is{Gam}] Chemical hardnesses of all species.

\item[\is{Kcn}] CN scaling factor of all species.

\item[\is{Rad}] Charge width of all species in Bohr.

\item[\is{Cutoff}] \modif{\modtype{length}} Cutoff distance when
  calculating electrostatics interactions under PBC.

\item[\is{EwaldParameter}] Sets the splitting parameter in the Ewald
  electrostatic summation for periodic calculations. This controls the fraction
  of the Ewald summation occurring in real and reciprocal space. Setting it to
  zero or negative activates an automatic determination of this parameter
  (default and recommended behaviour).

\item[\is{EwaldTolerance}] Sets the tolerance for Ewald summation in periodic
  systems.

\end{description}

\begin{verbatim}
ChargeModel = EEQ {
  EwaldParameter = 0.25165824
  EwaldTolerance = 1.0E-9
  Chi = Values {
    Ga = 1.15018618
    As = 1.36313743
  }
  Gam = Values {
    Ga = 8.299615E-2
    As = 0.19005278
  }
  Kcn = Values {
    Ga = -1.05627E-002
    As = 7.657769E-002
  }
  Rad = Values {
    Ga = 1.76901636
    As = 2.41244711
  }
}
\end{verbatim}

\subsubsection{DftD4 CoordinationNumber}
\label{sec:dftbp.CoordinationNumber}

The \kw{CoordinationNumber} determines how the local coordination environment
for its parent method is calculated. Currently four different counting functions
are available: \is{Erf\cb}, \is{Cov\cb}, \is{Gfn\cb}, and \is{Exp\cb}.
\is{Erf\cb} is the default coordination number for the EEQ charge model,
while \is{Cov\cb} is the default coordination number for DFTD4.

\begin{ptable}
  \kw{Electronegativities} & m & & \is{PaulingEN} & \\
  \kw{Radii} & m & & \is{CovalentRadiiD3} & \\
  \kw{Cutoff} & r & & 40 & \\
  \kw{CutCN} & r & & 0 / 8 & \\
\end{ptable}

\begin{description}

\item[\is{Radii}] Covalent radii of all species in Bohr.
  Default values taken are the DFTD3 covalent radii.\cite{grimme-jcp-132-154104}

\item[\is{Electronegativities}] Electronegativities of all species.
  Default values taken are Pauling ENs.

\item[\is{Cutoff}] \modif{\modtype{length}} Cutoff distance when
  evaluating counting function.

\item[\is{CutCN}] Maximum value for coordination number, coordination numbers
  higher than this value will be smoothly cut away. Deactivated for values
  smaller or equal to zero. Default depends on parent method.

\end{description}

\begin{verbatim}
CoordinationNumber = Cov {
  CutCN = 0
  Electronegativities = PaulingEN {}
  Radii = CovalentRadiiD3 {}
}
\end{verbatim}

\subsubsection{TS}
\label{sec:dftbp.DispTS}

Applies the Tkatchenko-Scheffler dispersion model~\cite{TkatchenkoPRL09,StohrJCP16}.

The following keywords can be used with this dispersion model:
\begin{ptable}
  \kw{EnergyAccuracy} & r & & 1e-7 & \\
  \kw{ForceAccuracy} & r & & 1e-6 & \\
  \kw{Damping} & r & & 20.0 & \\
  \kw{RangeSeparation} & & & 0.94 & \\
  \kw{ReferenceSet} & & & "TS" & \\
  \kw{ConvergentSCCOnly} & l & SCC = Yes & \is{Yes}\\
\end{ptable}
\begin{description}
\item[\is{EnergyAccuracy}] \modif{\modtype{energy}} Specifies the energy accuracy
  of the dispersion calculation.
\item[\is{ForceAccuracy}] \modif{\modtype{force}} Specifies the force accuracy
  of the dispersion calculation.
\item[\is{Damping}] Damping factor.
\item[\is{RangeSeparation}] Range separation parameter. The default value is the
  DFTB LDA-value. Note, that various DFTB-parametrisations may require different
  values.
\item[\is{ReferenceSet}] Reference values for the free atoms. Possible choices
  are \is{TS} (mostly used for molecules) and \is{TSsurf} (mostly used when
  calculating surfaces).
\item[\is{ConvergentSCCOnly}] If true, requires that the SCC cycle converges
  before evaluating TS dispersion {\it post-hoc}.
\end{description}

\subsubsection{MBD}
\label{sec:dftbp.DispMBD}

Applies the Tkatchenko many-body-dispersion model~\cite{AmbrosettiJCP14,StohrJCP16}.

The following keywords can be used with this dispersion model:
\begin{ptable}
  \kw{Beta} & r & & 0.83 & \\
  \kw{NOmegaGrid} & i & & 15 & \\
  \kw{KGrid} & 3r & &  & \\
  \kw{KGridShift} & & & 0.5 0.5 0.5 & \\
  \kw{VacuumAxis} & & & No No No & \\
  \kw{ReferenceSet} & & & "TS" & \\
  \kw{ConvergentSCCOnly} & l & SCC = Yes & \is{Yes}\\
\end{ptable}
\begin{description}
\item[\is{Beta}] Range separation parameter. Default value is optimised for pure
  DFT LDA functionals. The value may be different for various SK-sets.
\item[\is{NOmegaGrid}] \modif{\modtype{force}} Number of integration points.
\item[\is{KGrid}] Monkhorst-Pack-like k-point grid to integrate the Hamiltonian
  of the interacting dipoles over the Brillouin-zone. Note, that this grid is
  independent from the k-point grid used to integrate the atomic Hamiltonian.
\item[\is{KGrid}] Shift of the Monkhorst-Pack grid.
\item[\is{VacuumAxis}] Set \is{Yes} for each axis, where there is vacuum between
  the periodically repeated images of the central cell.
\item[\is{ReferenceSet}] Reference values for the free atoms. Possible choices
  are \is{TS} (mostly used for molecules) and \is{TSsurf} (mostly used when
  calculating surfaces).
\item[\is{ConvergentSCCOnly}] If true, requires that the SCC cycle converges
  before evaluating MBD dispersion {\itshape post-hoc}.
\end{description}

\subsection{DFTB3}
\label{sec:dftbp.DFTB3}
\index{DFTB3}

If you would like to use what is called ``DFTB3'' in some publication(s)
\cite{gauss-jctc-7-931}, this group of options include the relevant
modifications to the SCC Hamiltonian and energy. \emph{To enable the DFTB3
  model} you will need to set \is{{ThirdOrderFull} = Yes} and damp H--X the
interactions (see Section \ref{sec:dftbp.hcorr}).

\begin{description}

\item[\is{ThirdOrder}] If set to \is{Yes} the \textit{on-site} 3rd order
  correction \cite{yang-JPCA-111-10861} is switched on. This corrects the
  SCC-Hamiltonian with the derivatives of the Hubbard U parameters, which you
  have to specify for every element in \is{HubbardDerivs}. This correction only
  alters the on-site elements and is only maintained for backward
  compatibility. \emph{You should use the full version \is{ThirdOrderFull}
    instead.}

\item[\is{ThirdOrderFull}] If set to \is{Yes} the \textit{full} 3rd order
  correction \cite{gauss-jctc-7-931} is switched on. This corrects the
  SCC-Hamiltonian with the derivatives of the Hubbard U parameters, which you
  have to specify for every element in \is{HubbardDerivs}.

\item[\is{HubbardDerivs}] Derivatives of the Hubbard U for the 3rd order
  correction (on-site or full). For every element the appropriate parameter (in
  atomic units) must be specified. If you use shell resolved SCC (with full
  3rd order), you must specify a list of derivatives for every element, with one
  Hubbard U derivative for each shell of the given element.
\begin{verbatim}
Hamiltonian = DFTB {
  :
  ThirdOrder = Yes
  HubbardDerivs {
    O = -0.14
    H = -0.07
  }
  :
}
\end{verbatim}
\end{description}

\subsection{Implicit Solvation Model}
\label{sec:dftbp.Solvation}

\subsubsection{Generalized Born Model}
\label{sec:dftbp.GeneralizedBorn}

In generalized Born (GB) theory,\cite{onufriev2019} a molecule is considered as
continuous region with a dielectric constant~$\epsilon_\text{in}$ surrounded by
infinite solvent with a dielectric constant~$\epsilon_\text{out}$.
Charges~$q_\text{A}$ are located at the atomic sites~$\vec R_\text{A}$
and their interaction in the presence of a polarized solvent can be
expressed as the solvation energy
%
\begin{equation}
  \Delta G_\text{GB} =
  -\frac12 \left(\frac1{\epsilon_\text{in}}-\frac1{\epsilon_\text{out}}\right)
  \sum_{\text{A}=1}^N\sum_{\text{B}=1}^N\frac{q_\text{A}q_\text{B}}
  {\left(R^2_\text{AB} + a_\text{A}a_\text{B}\exp\left[-\frac{R^2_\text{AB}}
  {4a_\text{A}a_\text{B}}\right]\right)^{\frac12}}.
\end{equation}

where $a_\text{A/B}$ are the effective Born radii of the atoms A/B.
The GB model is added to the Hamiltonian as second order fluctuation
in the charge density, similar to the coulombic interactions.

The Born radii are evaluated by an Onufriev--Bashford--Case~(OBC) corrected
pairwise approximation to the molecular volume given as
%
\begin{equation}
  \frac1{a_\text{A}} = \frac1a_\text{scale}\left(
  \frac1{R^\text{cov}_\text{A} - R_\text{offset}}
  - \frac1{R^\text{cov}_\text{A}}
  \cdot \tanh\left[b_\text{OBC}\Psi_\text{A}
  - c_\text{OBC}\Psi_\text{A}^2 + d_\text{OBC}\Psi_\text{A}^3\right]
  \right)
\end{equation}

where $a_\text{scale}$ is a scaling factor for the Born radii,
$R_\text{offset}$ is a global shift parameter for the covalent radii and
$a/b/c_\text{OBC}$ are the coefficients for the volume polynomial in
the OBC correction to the Born radii.
$\Psi_\text{A}$ is the pairwise approximation to the volume integral
given by
%
\begin{equation}
  \Psi_\text{A} = \frac{R^\text{cov}_\text{A} - R_\text{offset}}2
  \sum_\text{B} \Omega(R_\text{AB}, R^\text{cov}_\text{A}, s_\text{B}R^\text{cov}_\text{B})
\end{equation}

with $\Omega$ being the pairwise function used to approximate the
volume integral, which is only dependent on the distance and the
covalent radii. Note that, the covalent radius of the second atom
is scaled by the element-specific descreening value~$s_\text{B}$
to compensate the systematic overestimation of the volume by this
approach.

To use the generalized Born model in the SCC procedure use the
\is{GeneralizedBorn\cb} method in the input to \kw{Solvation}.
The non-polar solvent area model can be combinded with the GB model
enabling to additionally correct for hydrogen bonding, the resulting
model is called GBSA.
The parameters for the GBSA model are currently available at
\url{https://github.com/grimme-lab/gbsa-parameters} and can be read in with
\kw{ParamFile} and will setup the complete \kw{GeneralizedBorn\cb} input.

Note that the GB(SA) model implemented is only available for finite
systems.

\begin{ptable}
  \kw{ParamFile} & s & & & \\
  \kw{Solvent} & m & not has \kw{ParamFile} & & \\
  \kw{FreeEnergyShift} & r & & & \\
  \kw{Temperature} & r & & 298.15 K & \\
  \kw{State} & s & & gsolv & \\
  \kw{Kernel} & s & & Still & \\
  \kw{BornScale} & r & & & \\
  \kw{BornOffset} & r & & & \\
  \kw{OBCCorrection} & 3r & & 1.00, 0.80, 4.85 & \\
  \kw{CM5} & m & & & \pref{sec:dftbp.CM5} \\
  \kw{Radii} & m & & \kw{vanDerWaalsRadiiD3} & \\
  \kw{Descreening} & m & & & \\
  \kw{Cutoff} & r & & 35 AA & \\
  \kw{SASA} & m & & & \pref{sec:dftbp.SolventAreaModel} \\
  \kw{HBondCorr} & l & has \kw{SASA} & & \\
  \kw{HBondStrength} & m & HBondCorr = Yes & & \\
  \kw{ALPB} & l & & No & \\
  \kw{Alpha} & r & ALPB = Yes & 0.571412 & \\
\end{ptable}

\begin{description}
  \item[\is{ParamFile}] Reads in a parameter file for GBSA, specifying this
    keyword automatically provides the \is{Solvent} information, and defaults
    values for \is{FreeEnergyShift}, \is{BornOffset}, \is{BornScale}, \is{SASA\cb}
    \is{Descreening}, \is{HBondCorr} and \is{HBondStrength}.
    Usually no other keywords need to be specified when \is{ParamFile} is present.

  \item[\is{Solvent}] Descriptors of the solvent, can be load from a database
    by providing the solvent name as string in the \is{FromName\cb} method or
    by specifying the necessary values with the \is{FromConstants\cb} method.
    \is{FromConstants\cb} requires the dielectric constant as real for \kw{Epsilon},
    the molecular mass in \kw{MolecularMass} \modif{\modtype{mass}} and the
    solvent density in \kw{Density} \modif{\modtype{massdensity}}.
    \kw{MolecularMass} and \kw{Density} only affect the calculation if
    ``reference'' is chosen as state of the solution.

  \item[\is{FreeEnergyShift}] \modif{\modtype{energy}} Shift for free energy
    calculations.

  \item[\is{Temperature}] Temperature for free energy calculations. Default is
    ambient temperature: 298.15 K.
    Only affects the calculation for if ``reference'' or ``mol1bar'' is chosen
    as state of the solution.

  \item[\is{State}] \modif{\modtype{energy}}
    Reference state of the solution for free energy calculations.
    The calculated state shift is added to the free energy shift.
    Takes ``gsolv'' (default), ``reference'' or ``mol1bar''.
    The reference state ``gsolv'' corresponds to 1\;l of ideal gas and 1\;l of
    liquid  solution,
    ``reference'' corresponds to 1\;bar of ideal gas and 1\;mol/L of liquid
    solution at infinite dilution,
    ``mol1bar'' corresponds to 1\;bar ideal gas and 1\;mol/L of liquid solution.

  \item[\is{Kernel}]
    Interaction kernel for the screened Coulomb operator.
    Possible options are the canonical ``Still'' kernel\cite{still1990} and the
    ``P16'' kernel.\cite{lange2012}

  \item[\is{BornScale}] Value for scaling of Born radii.

  \item[\is{BornOffset}] \modif{\modtype{length}}
    Offset value for Born radii calculation.

  \item[\is{OBCCorrection}] Parameters for Onufriev--Bashfold--Case volume
    polynomial to correct Born radii calculation.
    The default values 1.0, 0.8, 4.85 correspond to GB\textsuperscript{OBC}II,
    alternatively 0.8, 0.0, 2.91 can be used for GB\textsuperscript{OBC}I.\cite{onufriev2004}

  \item[\is{CM5}] Use the charge model 5 to correct the atomic partial charges
    before evaluating the Born energy.

  \item[\is{Radii}] Atomic radii for each element in Bohr, either takes
    \is{VanDerWaalsRadiiD3\cb} for DFT-D3 van-der-Waals radii (can be overwritten)
    or requires to provide \is{Values\cb} for all species.
    Both methods accept \modif{\modtype{length}} units.

  \item[\is{Descreening}] Descreening values for each species.
    Disabled by \is{Unity\cb} method or enabled by providing \is{Values\cb}
    for each species.

  \item[\is{Cutoff}] \modif{\modtype{length}}
    Real space cutoff for the calculation of the Born radii.

  \item[\is{HBondCorr}]
    Include an empirical hydrogen bond correction.
    Only available for GBSA models.

  \item[\is{HBondStrength}]
    Hydrogen bonding strength for each species used in the empirical hydrogen
    bond correction. To disable the correction for species not involved
    in hydrogen bonding, set the value to zero.

  \item[\is{ALPB}]
    Use analytical lineared Poisson Boltzmann (ALPB) model\cite{sigalov2006}
    instead of Generalized Born.
    The main difference is a molecular shape dependent contribution for charged
    systems.

  \item[\is{Alpha}]
    Alpha parameter in the ALPB model, should be kept to 0.571412 which was
    derived from first principles.\cite{sigalov2006}
\end{description}

Example for a GB model with CS2 as solvent:

\begin{verbatim}
Hamiltonian = DFTB {
  Solvation = GeneralizedBorn { # GFN2-xTB/GBSA(CS2)
    Solvent = fromName { "cs2" }
    FreeEnergyShift [kcal/mol] = 2.81072250
    Radii = vanDerWaalsRadiiD3 [Bohr] {}
    Descreening = Values {
       H = 0.93699367
       C = 0.83307834
       N = 1.02661619
       O = 0.96508008
    }
    BornScale = 1.40636177
    BornOffset [Bohr] = 1.653719965215E-03
    OBCCorrection = {1.00 0.80 4.85}
    Cutoff = 40
  }
}
\end{verbatim}

Example for a GBSA model for water

\begin{verbatim}
Hamiltonian = DFTB {
  Solvation = GeneralizedBorn { # GFN2-xTB/GBSA(water)
    Solvent = fromConstants {
      Epsilon = 80.2
      MolecularMass [amu] = 18.0
      Density [kg/l] = 1.0
    }
    FreeEnergyShift [kcal/mol] = 1.16556316
    BornScale = 1.55243817
    BornOffset = 2.462811043694508E-02
    Radii = vanDerWaalsRadiiD3 {}
    Descreening = Values {
      H = 0.71893869
      C = 0.74298311
      N = 0.90261230
      O = 0.75369019
    }
    SASA {
      ProbeRadius = 1.843075777670416
      Radii = vanDerWaalsRadiiD3 {}
      SurfaceTension = Values {
        H = -3.34983060E-01
        C = -7.47690650E-01
        N = -2.31291292E+00
        O =  9.17979110E-01
      }
    }
    HBondCorr = Yes
    HBondStrength = Values {
        H = -7.172800544988973E-02
        C = -2.548469535762511E-03
        N = -1.976849501504001E-02
        O = -8.462476828587280E-03
    }
  }
}
\end{verbatim}


\subsubsection{Charge Model 5}
\label{sec:dftbp.CM5}

The charge model 5 (CM5)\cite{marenich2012} can be used to correct partial charges by
%
\begin{equation}
   q^\text{CM5}_\text{A} = q_\text{A}
   + \sum_{B} D_\text{A--B} \exp[-\alpha(R_\text{AB} - R^\text{cov}_\text{A} - R^\text{cov}_\text{B}]
\end{equation}

The pairwise parameters $D_\text{A--B}$ are fixed to the original published ones,
while the exponent $\alpha$ and the covalent radii $R^\text{cov}_\text{A}$ default
to the published parameters but can be adjusted in the input.

\begin{ptable}
  \kw{alpha} & r & & 2.474 1/AA & \\
  \kw{Radii} & m & & \is{atomicRadii\cb} & \\
  \kw{Cutoff} & r & & 30.0 & \\
\end{ptable}

\begin{description}

  \item[\is{alpha}] \modif{\modtype{1/length}}
    Exponent of the CM5 correction.

  \item[\is{Radii}] Atomic covalent radii for each species, either takes
    \is{AtomicRadii\cb} for default atomic radii\cite{mantina2010}
    (can be overwritten) or requires to provide \is{Values\cb} for all species.
    Both methods accept \modif{\modtype{length}} units.

  \item[\is{Cutoff}] \modif{\modtype{length}}
    Real space cutoff for the calculation of the CM5 correction.
\end{description}

\begin{verbatim}
CM5 {
  Alpha = 1.30918451402600
  Radii = AtomicRad {}
}
\end{verbatim}


\subsubsection{Solvent area model}
\label{sec:dftbp.SolventAreaModel}

The non-polar solvation free energy can be estimated from the solvent accessible
surface area (SASA) by
%
\begin{equation}
  \Delta G_\text{non-polar} = \sum_\text{A} \gamma_\text{A} \sigma_\text{A}
\end{equation}
%
where $\gamma_\text{A}$ is the surface tension and $\sigma_\text{A}$
is the accessible surface area of each atom.
To calculate the latter, the molecule is assumed as a convolution of spheres
which is probed by a probe sphere rolled around the surface.
Here a smooth numerical integration approach is employed.\cite{im2003}

To use the non-polar surface area model in an calculation use the
\is{SASA\cb} method in the input to \kw{Solvation}.
This model is currently only available for finite systems.

\begin{ptable}
  \kw{ProbeRadius} & r & & & \\
  \kw{Smoothing} & r & & 0.3 AA & \\
  \kw{Offset} & r & & 2.0 AA & \\
  \kw{Tolerance} & r & & 1.0e-6 & \\
  \kw{AngularGrid} & i & & 230 & \\
  \kw{Radii} & m & & \kw{vanDerWaalsRadiiD3} & \\
  \kw{SurfaceTension} & m & & & \\
\end{ptable}

\begin{description}
  \item[\is{ProbeRadius}] \modif{\modtype{length}}
    Radius of the probe sphere used to determine the accessible surface area.

  \item[\is{Smoothing}] \modif{\modtype{length}}
    Smoothing parameter for numerical integration.

  \item[\is{Offset}] \modif{\modtype{length}}
    This offset value is added on the realspace cutoff radius for the
    neighbourlist generation.
    The realspace cutoff is determined automatically from the probe radius,
    the largest atomic radius and the smoothing parameter.

  \item[\is{Tolerance}]
    Minimal value of surface area contribution of a grid point to be accounted
    for as SASA.

  \item[\is{AngularGrid}]
    Size of the angular Lebedev--Laikov integration grid.\cite{lebedev1999}
    The grid size mainly determines the computational cost of evaluating the
    accessible surface area, too small grid sizes can lead to significant errors
    due to missing rotational invariance.
    A safe choice should be 230 grid points per atom.
    Possible values are
    6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, \textit{230}, 266, 302, 350,
    434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890,
    4334, 4802, 5294, 5810.

  \item[\is{Radii}] Atomic radii for each element, either takes
    \is{VanDerWaalsRadiiD3\cb} for DFT-D3 van-der-Waals radii (can be overwritten)
    or requires to provide \is{Values\cb} for all species.
    Both methods accept \modif{\modtype{length}} units.

  \item[\is{SurfaceTension}]
    Surface tension parameter for each species in dyn/cm.

\end{description}

\begin{verbatim}
Hamiltonian = DFTB {
  Solvation = SASA { # GFN1-xTB/GBSA(Toluene)
    ProbeRadius [AA] = 1.59772343
    Smoothing [AA] = 0.3
    Offset [AA] = 2
    AngularGrid = 230
    Radii = vanDerWaalsRadiiD3 {}
    SurfaceTension = Values {
      H = -1.52312760
      C = -2.92375089
      O = 0.79482640
    }
  }
}
\end{verbatim}


\subsection{Halogen corrections}
\label{sec:dftbp.xcorr}

The \kw{HalogenXCorr} keyword includes the halogen correction of
Ref.~\cite{kubillus-jctc-11-332}. This is fitted for the DFTB3-D3 model and the
{\tt 3ob-3-1} parameter set. The correction is only relevant for systems
including interactions between \{O,N\}--\{Cl,Br,I\} pairs of atoms.


\subsection{Hydrogen corrections}
\label{sec:dftbp.hcorr}

There are currently two available methods to correct hydrogen interactions
(mainly hydrogen bonds) in the \is{HCorrection} environment:

\subsubsection{Damping}
\index{sec:Damp X-H}

The \is{Damping} method modifies the short range contribution to the SCC
interaction between atoms $A$ and $B$ with the damping factor
\begin{equation*}
  e^{-\left(\frac{U_{Al} + U_{Bl}}{2}\right)^\zeta r_{AB}^2}
\end{equation*}
provided that at least one of the two atoms is
hydrogen~\cite{gauss-jctc-7-931,yang-JPCA-111-10861}. ($U_{Al}$ and $U_{Bl}$ are
the Hubbard U values of the two atoms for the $l$-shell, $r_{AB}$ is the
distance between the atoms.) An atom is considered to be a hydrogen-like atom,
if its mass (stored in the appropriate homonuclear SK-file) is less than 3.5
amu.  The \is{Exponent} keyword in this environment sets the parameter $\zeta$
for the short range damping:
\begin{verbatim}
HCorrection = Damping {
  Exponent = 4.05
}
\end{verbatim}
Table 2 of reference~\cite{gauss-jctc-7-931} gives suggested values of the
exponent for different DFTB2 and DFTB3 models applied to light atoms bonded to
hydrogen.

\subsubsection{DFTB3-D3H5}
\index{sec:DFTB3-D3H5}

DFTB3-D3H5~\cite{rezac-jctc-13-2017} is a variant of DFTB3 with additional
corrections for non-covalent interactions (dispersion and hydrogen bonds).  It
consists of a third-order DFTB calculation using the 3OB parameter set, but
where the gamma-function damping (\is{Damping} method above) is replaced by the
H5 correction and an additional D3 dispersion correction in included. This
method also includes a repulsive term which is added to prevent unphysically
close approach of pairs of hydrogen atoms~\cite{rezac-jctc-8-2012}.

Setting the \is{HCorrection} environment to \iscb{H5} activates this
correction for hydrogen bonds~\cite{rezac-jctc-13-2017}. If no additional
parameters are provided in the input, suitable values for H-\{O,N,S\} systems
are used (the correction was developed for the DFTB3/3OB model and parameters).
\begin{verbatim}
HCorrection = H5 {}
\end{verbatim}

\textbf{Note:} It was found that DFTB3 overestimates the strength of H-bonds
involving the terminal nitrogen of an azide group, and the published results in
Ref.~\cite{rezac-jctc-13-2017} were obtained with the H5 correction switched off
for these specific atoms. To reproduce this behavior in a system containing
nitrogen in several environments, a new atom type with a different name but the
same DFTB parameters can be used for specific N atoms to which the correction
should not be applied.

If you want to specify the parameters manually, \is{H5} accepts following
options, corresponding to terms in Ref.~\cite{rezac-jctc-13-2017}:
\begin{ptable}
  \kw{RScaling} & r & & 0.714 & \\
  \kw{WScaling} & r & & 0.25 & \\
  \kw{H5Scaling} & m & &  & \\
\end{ptable}
\begin{description}
\item[\is{RScaling}] Global scaling factor, $s_r$, when calculating the position
  of the correcting gaussian functions:
  \begin{equation*}
    r_0 = s_r \left(r_\text{vdW}(X) + r_\text{vdW}(H)\right) \text.
  \end{equation*}

\item[\is{WScaling}] Global scaling factor, $s_W$, when calculating the width of
  the correcting gaussian functions. The full-width at at half-maximum of the
  gaussian, $w$, is normalised to be 1 for a unit value of \is{WScaling}:
  \begin{equation*}
   w = \frac{s_w \left(r_\text{vdW}(X) + r_\text{vdW}(H)\right)}{2\sqrt{2 \ln
       2}} \text.
  \end{equation*}

\item[\is{H5Scaling}] Atom type specific scaling pre-factor, $k_{X\text{H}}$, of
  the correcting gaussian functions when calculating the SCC-interaction:
  \begin{equation*}
    \gamma_{X\text{H}}^{\text{H5}} = \gamma_{X\text{H}} \left( 1 + k_{X\text{H}}
    \exp\left(-\frac{\left(r_{X\text{H}} - r_0\right)^2}{2w^2}\right)\right)
    \text.
  \end{equation*}
  You will have to specify one value for each of the chemical species you would
  like to correct (see the example below). Explicitly setting a negative value
  (e.g. \is{-1.0}) for a given atom type switches off the correction for
  hydrogen bonds involving that type of atom. In the special cases of N, O or S,
  if you do not specify a value (and do not disable the contribution by using
  \is{-1.0}), the default value from the reference paper will be
  used~\cite{rezac-jctc-13-2017}. For any other omitted atom types, the code
  defaults to a choice of \is{-1.0} (no correction).
\end{description}

\begin{verbatim}
Hamiltonian = DFTB {
  :
  HCorrection = H5 {
    RScaling = 0.714
    WScaling = 0.25
    H5Scaling {
      O = 0.06
      N = 0.18
      S = 0.21
    }
  }
  :
}
\end{verbatim}

\textbf{Note:} \label{page:dftbp.H5} The van der Waals radii ($r_\text{vdW}$) of
atoms are also required. \dftbp{} stores these for most of the periodic table,
but for cases that are not available their contribution to this correction are
neglected.

For a DFTB3-D3H5 calculation, a specific parametrization of the D3 dispersion
has to be used. In addition to setting up appropriate values of the D3
parameters, as discussed in Ref.~\cite{rezac-jctc-13-2017}, the
hydrogen--hydrogen repulsion of Ref.~\cite{rezac-jctc-8-2012} has to also be
activated. The complete input is:
\begin{verbatim}
Hamiltonian = DFTB {
  :
  Dispersion = DftD3 {
    Damping = ZeroDamping {
      sr6 = 1.25
      alpha6 = 29.61
    }
    s6 = 1.0
    s8 = 0.49
    HHRepulsion = Yes
  }
  :
}
\end{verbatim}

\subsection{RangeSeparated}
\label{sec:dftbp.RangeSep}

The \kw{RangeSeparated} keyword specifies the use of a range separated hybrid
functional. Currently, only the long-range corrected hybrid functional (\is{LC})
\cite{niehaus-PSSB-249-237,lutsker-JCP-143-184107} is implemented. There, the
electrostatic interaction is split up into long and short ranged components
according to
\begin{equation*}
\frac{1}{r}=\frac{1-e^{-\omega r}}{r}+\frac{e^{-\omega r}}{r},
\end{equation*}
with the range-separation parameter $\omega$, which is set in the Slater-Koster
files. The option should only be used with corresponding parameter sets created
for use with long-range corrections.

The \kw{RangeSeparated} keyword expects either \is{None} (default -- no use of
range-separated hybrid functional) or the \is{LC}\cb\ block as value. The latter
has the following option:

\begin{ptable}
  \kw{Screening} & m & & Thresholded \cb & \\
\end{ptable}


\begin{description}
\item[\is{Screening}] Choice of the screening method. The following choices are
  possible:

  \begin{description}
  \item[\is{Thresholded \{\}}] Screening according to estimated magnitude of
    terms. This is faster than the \is{NeighbourBased} method below, but does
    not support all of the cases (restarting and spin polarisation).

    \begin{ptable}
      \kw{Threshold} & r & & 1e-6 & \\
      \kw{CutoffReduction} & r & & 0.0 & \\
    \end{ptable}

    \begin{description}
    \item[\is{Threshold}] Threshold, below which elements are considered to be
      zero.
    \item[\is{CutoffReduction}]\modif{\modtype{length}} Reduces the spatial
      cutoff, beyond which the overlap between atoms is considered to be zero.
      This can be used as an additional tweak to speed up the LC-calculation,
      but make sure first, that your results do not change
      considerably. Default: 0.0 -- no reduction, using the cutoff from the
      SK-files.
    \end{description}

  \item[\is{NeighbourBased}] Uses a purely neighbour-list based algorithm. This
    algorithm is usually considerably slower than the \is{Thresholded}.

    \begin{ptable}
      \kw{CutoffReduction} & r & & 0.0 & \\
    \end{ptable}

    \begin{description}
    \item[\is{CutoffReduction}]\modif{\modtype{length}} See description in the
      \is{Thresholded} block.
    \end{description}

  \item[\is{MatrixBased}] Uses a matrix-matrix multiplication based algorithm. This
    can be faster than other two algorithms.

  \end{description}
\end{description}

Example for thresholded screening with customised threshold value.
\begin{verbatim}
RangeSeparated = LC {
  Screening = Thresholded {
    Threshold = 1e-5
  }
}
\end{verbatim}

Example for neighbour list based screening with customised cutoff reduction:
\begin{verbatim}
RangeSeparated = LC {
  Screening = NeighbourBased {
    CutoffReduction [AA] = 2.0
  }
}
\end{verbatim}

Example for matrix-matrix multiplication based method:
\begin{verbatim}
RangeSeparated = LC {
  Screening = MatrixBased {}
}
\end{verbatim}


\subsection{On site corrections}
\label{sec:dftbp.Onsites}

This block enables corrections for on-site matrix elements which improve the
description of multi-centre integrals~\cite{garcia14Thesis} leading to, for
example, improved hydrogen-bond energies~\cite{dominguez15}.

For each chemical species, the spin-same-spin and spin-different-spin constants
should be specified for all combinations of atomic shells. {\bf note:} the
matrix of constants is symmetric and the purely s-with-s entries are zero (the
code ignores their value due to symmetry).

Example:
\begin{verbatim}
OnSiteCorrection= {
    # same spin oxygen
    Ouu = {0.00000  0.08672
           0.08672 -0.00523}
    # hetero-spin oxygen
    Oud = {0.00000  0.14969
           0.14969  0.03834}
    # H all zero
    Huu = {0}
    Hud = {0}
}
\end{verbatim}

Some on-site constants are given in appendix~\ref{app:onsiteconst}.

\subsection{Differentiation}
\label{sec:dftbp.Differentiation}

Calculations of forces currently require the numerical derivatives of the
overlap and non-self-consistent Hamiltonian. This environment controls how these
derivatives are evaluated.

\textbf{Note:} In earlier DFTB+ versions (up to version 1.2), differentiation
was done using finite difference derivatives with a step size of 0.01 atomic
units. If you want to reproduce old results, choose the \is{FiniteDiff} method
and set the step size explicitly to this value.

\subsubsection{FiniteDiff\cb}
\label{sec:dftbp.FiniteDiff}

Finite difference derivatives with a specified step size
\begin{ptable}
  \kw{Delta} & r & & epsilon$^{\sfrac{1}{4}}$ & \\
\end{ptable}
\begin{description}
\item[\is{Delta}]\modif{\modtype{length}} Step size
\end{description}


\subsubsection{Richardson\cb}
\label{sec:dftbp.Richardson}

Extrapolation of finite difference via Richardson's deferred approach to the
limit (in principle the most accurate of the currently available choices).

\subsection{ForceEvaluation}
\label{sec:dftbp.ForceEvaluation}

Chooses the method for evaluating the electronic contribution to the
forces.

\begin{description}
\item[\is{'traditional'}] Uses the ``traditional'' DFTB-force expression, given
  for example, in Ref.~\cite{elstner-prb-58-7260}.
\item[\is{'dynamics'}] Force expression from
  Ref.~\cite{aradi-jctc-11-3357}. This choice should be used if forces are being
  calculated with non-converged charges (e.g.\ when doing XLBOMD
  dynamics). \textbf{Note:} this force expression is only compatible with the
  Fermi filling (see keyword \is{Filling}, p.~\pref{sec:dftbp.Filling}.)
\item [\is{'dynamicsT0'}] Simplified dynamic force expression valid for
  electronic temperature \mbox{$T=0$~K} \cite{aradi-jctc-11-3357}.  This choice
  should be used if forces are calculated with non-converged charges and the
  electronic temperature is zero (e.g.\ when doing XLBOMD dynamics at
  \mbox{$T=0$~K}).
\end{description}

\textbf{Note:} that XLBOMD calculations (Section \ref{sec:dftbp.xlbomd}) are not
able to use the \is{'traditional'} forces.

Example:
\begin{verbatim}
ForceEvaluation = 'dynamics'
\end{verbatim}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Options}
\label{sec:dftbp.Options}

This block collects some global options for the run.
\begin{ptable}
  \kw{WriteAutotestTag} &l & & No & \\
  \kw{WriteDetailedXML} & l & & No & \\
  \kw{WriteResultsTag} & l & & No & \\
  \kw{WriteDetailedOut} & l & & Yes & \\
  \kwl{RestartFrequency}{kw:dftbp.RestartFrequency} & i &Driver = \cb, SCC = Yes & 20 & \\
  \kw{RandomSeed} & i & & 0 & \\
  \kw{MinimiseMemoryUsage} & l & & No & \\
  \kw{TimingVerbosity} & i & & 0 & \\
  \kw{ShowFoldedCoords} & l & Periodic = Yes & No & \\
  \kwl{WriteHS}{kw:dftbp.WriteHS} & l & & No & \\
  \kwl{WriteRealHS}{kw:dftbp.WriteRealHS} & l & & No & \\
  \kwl{ReadChargesAsText}{kw:dftbp.ReadChargesAsText} & l &
  ReadInitialCharges = Yes & No & \\
  \kwl{WriteChargesAsText}{kw:dftbp.WriteChargesAsText} & l & & No & \\
  \kw{SkipChargeTest} & l & ReadInitialCharges = Yes & No & \\
\end{ptable}
\begin{description}

\item[\is{WriteAutotestTag}] Turns the creation of the
  \verb|autotest.tag| file on and off. (This file can get quite big
  and is only needed for the autotesting framework.)

\item[\is{WriteDetailedXML}] Turns the creation of the
  \verb|detailed.xml| file on and off. (The \verb|detailed.xml| file
  is needed among others by the \verb|waveplot| utility for
  visualising molecular orbitals.)

\item[\is{WriteResultsTag}] Turns the creation of the \verb|results.tag| file on
  and off. (That file is used by several utilities processing the results of
  \dftbp.) For a description of the file format see
  p.~\pref{sec:dftbp.results}.

\item[\is{WriteDetailedOut}] Controls the creation of the file
  \verb|detailed.out| (see p.~\pref{sec:dftbp.detailedout}). Since this contains
  the detailed information about the last step of your run, you shouldn't turn
  it off without good reasons.

\item[\is{RestartFrequency}] Specifies the interval at which charge
  restart information should be written to disc for static SCC
  calculations. Setting it to \is{0} prevents the storage of restart
  information. If running an MD calculation, see also
  section~\ref{kw:dftbp.MDRestartFrequency} regarding \is{MDRestartFrequency}.

\item[\is{RandomSeed}] Sets the seed for the random number
  generator. The value \is{0} causes random initialisation. (This
  value can be used to reproduce earlier MD calculations by setting
  the initial seed to the same value.)

\item[\is{MinimiseMemoryUsage}] Tries to minimise memory usage by
  storing various matrices on disc instead of keeping them in memory.
  Set it to \is{Yes} to reduce the memory requirement for calculations
  with many k-points or spin polarisation. Note: Currently this option has no
  effect and you will get a warning if setting it to be \is{Yes}.

\item[\is{TimingVerbosity}] Level of information regarding CPU and wall clock
  timings of sections of the code, higher values becoming more verbose. Setting
  this parameter to 0 or below supresses any information being printed
  (default). Setting it to -1 includes all measured timings.

\item[\is{ShowFoldedCoords}] Print coordinates folded back into the
  central cell, so if an atom moves outside the central cell it will
  reappear on the opposite side. The default behaviour is to use
  unfolded coordinates in the output. (Please note, that this option
  only influences how the coordinates are printed and written, it does
  not change the way, periodic systems are treated internally.)

\item[\is{WriteHS}] Instructs the program to build the square
  Hamiltonian and overlap matrices and write them to files. The output
  files are \verb|hamsqrN.dat| and \verb|oversqr.dat|, where \verb|N|
  enumerates the spin channels. For a detailed description of the file
  format see p.~\pref{sec:dftbp.hamsqr}.

  \textbf{Note:} If either of the options \is{WriteHS} or \is{WriteRealHS} are
  set to \is{Yes}, the program only builds the matrices, writes them to disc and
  then stops immediately. No diagonalisation, no SCC-cycles or geometry
  optimisation steps are carried out. You can use the \is{ReadInitialCharges}
  option to build the Hamiltonian with a previously converged charge
  distribution.

\item[\is{WriteRealHS}] Instructs the program to build the real space
  (sparse) Hamiltonian and overlap matrices and write them to
  files. The output files are \verb|hamreal.dat| and
  \verb|overreal.dat|. For a detailed description of the file format
  see p.~\pref{sec:dftbp.hamsqr}.

  \textbf{Note:} If either of the options \is{WriteHS} or \is{WriteRealHS} are
  set to \is{Yes}, the program only builds the matrices, writes them to disc and
  then stops immediately. No diagonalisation, no SCC-cycles or geometry
  optimisation steps are carried out. You can use the \is{ReadInitialCharges}
  option to build the Hamiltonian with a previously converged charge
  distribution.

\item[\is{ReadChargesAsText}] If \is{No}, the program expects the file
  \verb|charges.bin| to contain starting charges stored in binary. If \is{Yes},
  then \verb|charges.dat| should contain a text file of this data. See section
  \ref{sec:charges.bin}.

\item[\is{WriteChargesAsText}] If \is{No}, the program stores charges in the
  binary file \verb|charges.bin|, while if \is{Yes} then \verb|charges.dat|
  contains text of this data. See section \ref{sec:charges.bin}.

\item[\is{SkipChargeTest}] If \is{Yes}, testing of whether the charges read from
  file match the total charge (and magnetisation) specified in the {\dftbp}
  input (if relevant) is performed. Skipping this test (setting to \is{No}) may
  be useful if restarting from a charges generated for a similar system with
  slightly different total charge or magnetisation. Similarly, in the event of
  serious instabilities in the SCC cycle, the generated charge restart file may
  fall outside of the check-sum tolerances, hence this option allows a
  re-start. Finally, in the case of user edited \verb|charges.dat| file (see
  section~\ref{sec:charges.bin}), the check-sum this option removed the
  requirement that the checksum values in the file match the charges.

\end{description}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Analysis}
\label{sec:dftbp.Analysis}

This block collects some options to analyse the results of the
calculation and/or calculate properties.
\begin{ptable}
  \kw{AtomResolvedEnergies} & l & & No & \\
  \kw{MullikenAnalysis} & l & & Yes & \\
  \kw{CM5} & m & MullikenAnalysis = Yes & & \pref{sec:dftbp.CM5} \\
  \kw{WriteNetCharges} & l & MullikenAnalysis = Yes & No & \\
  \kw{ProjectStates} & m & & \cb & \\
  \kw{Localise} & m & & \cb & \\
  \kwl{WriteEigenvectors}{kw:dftbp.WriteEigenvectors} & l & & No & \\
  \kw{EigenvectorsAsText} & l & WriteEigenvectors = Yes & No & \\
  \kw{WriteBandOut} & l & & Yes & \\
  \kw{CalculateForces} & l & & No & \\
  \kw{ElectrostaticPotential} & m & SCC = Yes & \cb & \pref{sec:dftbp.esp} \\
\end{ptable}

\begin{description}

\item[\is{AtomResolvedEnergies}] Specifies whether the contribution of the
  individual atoms to the total energies should be displayed or not.

\item[\is{MullikenAnalysis}] If \is{Yes}, the results of a Mulliken analysis of
  the system is given.

\item[\is{WriteNetCharges}] If \is{Yes}, the net charges (component of the
  charge associated only with the on-site part of the atomic orbitals, not the
  off-site hybridisation).

\item[\is{CM5}] If present the charge model 5 (CM5)\cite{marenich2012} corrected
  atomic partial charges will be written.

\subsubsection{ProjectStates} \kw{ProjectStates}
\label{sec:dftbp.ProjectStates} evaluates the Mulliken
projection of electronic states onto specific regions of the system being
modelled (partial density of states -- PDOS). The format of the projected data
files is similar to \verb|band.out|, but the second column is the fraction of
the state within that region, instead of its occupation number (for
non-collinear and spin-orbit calculations, three additional columns for the
magnetisation of the state are also given).

Each region for projection is specified within a \iscb{Region} block,
with the following options
\begin{ptable}
  \kw{Atoms} & (i|s)+ &  & - & \\
  \kw{ShellResolved}  & l & & \is{No} & \\
  \kw{OrbitalResolved} & l & & \is{No} & \\
  \kw{Label} & s &  & "region{\it i}" & \\
\end{ptable}
\begin{description}
\item[\is{ShellResolved}] Project onto separate atomic shells of the
  region. These are taken in order of increasing shell number of the
  atoms. \is{ShellResolved = Yes} is only allowed, if all the selected atoms
  are of the same type.
\item[\is{OrbitalResolved}] Project onto separate atomic orbitals of the
  region. These are taken in order of increasing shell number of the atoms. As
  with \is{ShellResolved}, this only allowed, if all the selected atoms are of
  the same type.
\item[\is{Atoms}] Specification of the atoms over which to make the
  projection.. Atoms are specified in the same way as \is{MovedAtoms}
  in section \ref{sec:dftbp.SteepestDescent}.)
\item[\is{Label}] Prefix of the label for the resulting file of data for this region. The default is
  ``region\textit{i}.out'' where {\it i} is the number of the region in the input. In the case that
  \is{ShellResolved} = \is{Yes}, the shell index is appended, so that files with names
  ``Label.\textit{j}.out'' are written. For \is{OrbitalResolved = Yes}, the shell and then
  $m$-value is appended, so that files with names ``Label.\textit{j}.\textit{m}.out'' are written.
\end{description}
Examples:
\begin{verbatim}
  ProjectStates = {
    Region = {              # first region
      Atoms = 23:25 27      # atoms 23, 24, 25 and 27
    }
    Region = {
      Atoms = N             # All nitrogen atoms
      ShellResolved = Yes   # s and p shells separated instead of atomic PDOS
      Label = "N"           # files N.1.out and N.2.out for s and p states
    }
  }
\end{verbatim}

\subsubsection{Localise}
\label{sec:dftbp.Localise} Convert the single particle states of the calculation
to localised orbitals via a unitary transformation. Localised orbitals span the
same states as the occupied orbitals, so are equivalent to the usual valence
band states, but are more localised in space.  Currently only \kw{PipekMezey}
localisation is supported (but not for non-collinear or spin-orbit calculations).

Pipek-Mezey~\cite{pipek-JCP-90-4916} localisation transforms the occupied
orbitals such that the square of the Mulliken charges for each orbital is
maximised. The resulting localised states are output as
\verb|localOrbs.out|\index{localOrbs.out} and
\verb|localOrbs.bin|\index{localOrbs.bin} following the format given in
appendix~\ref{sec:dftbp.eigenvec} for \verb|eigenvec.out| and \verb|eigenvec.bin|.

\begin{ptable}
  \kw{Tolerance} & r &  & 1E-4 & \\
  \kw{MaxIterations} & i &  & 100 & \\
\end{ptable}
\begin{description}
\item[\is{Tolerance}] Cut off for rotations in the localisation process.
\item[\is{MaxIterations}] Maximum number of total sweeps to perform.
\end{description}

For systems with non-gamma-point $k$-points, no further options are available.
\begin{verbatim}
Analysis = {
  Localise = {
    PipekMezey = {
       # These are the default options, which are also set if the bracket is left empty.
       Tolerance = 1.0E-4
       MaxIterations = 100
    }
  }
}
\end{verbatim}

For molecular and gamma point periodic calculations there are two
implementations available, \is{Dense} = \is{Yes} will use the $O(n^4)$ scaling
conventional algorithm, while \is{Dense} = \is{No}, uses the default sparse
method which {\em may} have better scaling properties.

\begin{ptable}
  \kw{Dense} & l & & \is{No} & \\
  \kw{SparseTolerances} & r+ & Dense = No & 1E-1 1E-2 1E-6 1E-12 & \\
\end{ptable}
\begin{description}
\item[\is{Dense}] Selects the conventional method (\is{Yes}) using Jacobi sweeps
  over all orbital pairs or (\is{No}) uses the default sparse method.
\item[\is{SparseTolerances}] The sparse method introduces support regions during
  evaluation to increase performance, and these requires a set of tolerances to
  determine the regions to be used (these are listed in decreasing order, i.e.,
  with tighter tolerances as the localisation proceeds).
\end{description}

\item[\is{WriteEigenvectors}] Specifies, if eigenvectors should be printed in
  \verb|eigenvec.bin|. For a description of the file format see
  p.~\pref{sec:dftbp.eigenvec}.

\item[\is{EigenvectorsAsText}] If eigenvectors are being written, specifies if a
  text version of the data should be printed in \verb|eigenvec.out|. For a
  description of the file format see p.~\pref{sec:dftbp.eigenvec}.

\item[\is{WriteBandOut}]  Controls the creation of the file \verb|band.out|
  which contains the band structure in a more or less human friendly format.

\item[\is{CalculateForces}] If \is{Yes}, forces are reported, even if not needed
  for the actual calculation (e.g.\ static geometry calculation).

\end{description}

\subsubsection{ElectrostaticPotential}
\label{sec:dftbp.esp}

Evaluates the electrostatic potential at specified points in space for SCC
calculations. This data is accumulated in a specified text file.

\begin{ptable}
  \kw{OutputFile} & s  &  & "ESP.dat" & \\
  \kw{AppendFile} & l  & MD or geometry optimisation & No & \\
  \kw{Softening}  & r     & & 1E-6 & \\
  \kw{Points}     & (3r)+ & \is{Grid} not set & \cb  & \\
  \kw{Grid}       & m     & \is{Points} not set & \cb  & \\
\end{ptable}
\begin{description}
\item[\is{OutputFile}] Text file to store the potential. If external electric
  fields are present, an additional column gives their values. See
  p.~\pref{sec:dftbp.espfile} for a description of the file.
\item[\is{AppendFile}] If running calculations with multiple geometries, should
  the \is{OutputFile} be appended or only contain the last potential
  information?
\item[\is{Softening}] \modif{\modtype{length}} Modifies the plotted potential to
  remove the $r=0$ divergence of $\sfrac{1}{r}$, by setting $\epsilon$ and
  instead plotting $\sfrac{1}{\sqrt{r^2 +\epsilon^2}}$. Internal potential
  calculations are unaffected, only the exported data.
\item[\is{Points}] \modif{\modtype{length}} List of cartesian points at which to
  evaluate the electrostatic field. In the case Periodic = Yes, the modifier "F"
  may instead be used to specify the points as fractions of the lattice vectors.
\item[\is{Grid}] \modif{\modtype{length}} Specification of a regular 1, 2 or 3
  dimensional grid of points. In the case Periodic = Yes, the modifier "F" may
  instead be used to specify the points as fractions of the lattice vectors.
  \begin{ptable}
    \kw{GridPoints} & 3i  &  & & \\
    \kw{Origin} & 3r  &  & & \\
    \kw{Spacing} & 3r  &  & & \\
    \kw{Directions} & 9r  & Modifier not F & 1 0 0\ \  0 1 0\ \  0 0 1& \\
  \end{ptable}
  \begin{description}
    \item[\is{Spacing}] Separation between points in each direction. This
      inherits the modifier for \is{Grid}.
  \item[\is{Origin}] Location of first point in the grid. This inherits the
    modifier for \is{Grid}.
  \item[\is{GridPoints}] Number of points in each of the three direction of the
    grid (a value of 1 places all points at the \is{Origin} of that direction).
  \item[\is{Directions}] Set of 3 cartesian vectors along which the grid will
    become aligned. This can rotate, skew, {\it etc.}\ the grid. The vectors are
    internally normalised, but must be independent.
  \end{description}
\end{description}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Excited State
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ExcitedState}
\label{sec:dftbp.ExcitedState}

This block collects some options to calculate excited states.
\begin{ptable}
  \kw{Casida} & p & SCC = Yes & \cb & \\
  \kw{PP-RPA} & p & SCC = Yes & \cb & \\
\end{ptable}

\subsection{Casida}
\label{sec:dftbp.Casida}

This tag contains the specifications for a time-dependent DFTB calculation,
based on linear response theory~\cite{niehaus-prb-63-085108}.

\textbf{Note:} the \dftbp{} binary must be compiled with linear response
calculations enabled to make use of these features (the
ARPACK~\cite{Lehoucq97arpackusers} library or ARPACK-ng~\cite{ARPACK-ng} is
required).

The calculation of vertical excitation energies and the corresponding oscillator
strengths as well as excited state geometry optimisation can be performed with
these options, details of the resulting output files are given in appendix
\ref{sec:tddftb_lr}. Linear response theory is currently implemented only for
the SCC-DFTB level of theory and molecular systems.\footnote{Excitation energies
  can also be calculated for gamma point periodic systems, but will be incorrect
  for delocalised excitations or for charge transfer-type excited states.}
Excitations can be calculated for fractional occupations and collinear
spin-polarisation. Forces (and hence geometry optimisation or MD) are
available for collinear spin-polarised systems but not for fractional occupations. The
specifications for this block have the following properties:

\textbf{Note:} Excited state calculations with the \kw{RangeSeparated}
functionals (Section \ref{sec:dftbp.RangeSep}) can be performed for spin free
molecular (non-periodic) singlet systems that do not have fractional fillings.

  \begin{ptable}
    \kw{NrOfExcitations}        & i & & - & \\
    \kw{StateOfInterest}        & i & & 0 & \\
    \kw{Symmetry}               & s & SpinPolarisation $=$ \cb~ & - & \\
    \kw{EnergyWindow}           & r & & FORTRAN HUGE() & \\
    \kw{OscillatorWindow}       & r & & -1 & \\
    \kw{WriteTransitions}       & l & & No & \\
    \kw{WriteTransitionCharges} & l & \parbox{0.3\textwidth}{\raggedright
      \is{RangeSeparated} and \is{StateOfInterest} non-zero} & No & \\
    \kw{WriteSPTransitions}     & l & & No & \\
    \kw{WriteMulliken}          & l & & No & \\
    \kw{WriteCoefficients}      & l & & No & \\
    \kw{WriteEigenvectors}      & l & & No & \\
    \kw{WriteDensityMatrix}     & l & & No & \\
    \kw{TotalStateCoeffs}       & l & WriteCoefficients = Yes & No & \\
    \kw{WriteXplusY}            & l & & No & \\
    \kw{WriteTransitionDipole}  & l & & No & \\
    \kw{WriteStatusArnoldi}     & l & & No & \\
    \kw{TestArnoldi}            & l & & No & \\
    \kw{ExcitedStateForces}     & l & \is{CalculateForces} = Yes & Yes & \\
    \kw{CacheCharges}           & l & & Yes & \\
  \end{ptable}

  \begin{description}

  \item[\is{NrOfExcitations}] Specifies the number of vertical excitation
    energies to be computed for every symmetry (singlet or triplet). It is
    recommended that a value slightly greater than the actual number of the
    states of interest is specified (the eigenvalue solver may not converge to
    the right roots otherwise).

  \item[\is{StateOfInterest}] Specifies the target excited state or states that
    should be calculated. These are numbered from the first (lowest) excited
    state as 1, and so on. If the absorption spectrum at a given geometry is
    required (i.e., a single-point calculation), this parameter should be set to
    zero (default) and the \is{Driver} section (\ref{sec:dftbp.Driver}) should
    be left empty (forces will not be available). A value less than 0 requests
    that the state with the largest dipole transition moment be found (again a
    single-point calculation).

  \item[\is{Symmetry}] Specifies the spin symmetry of the excited states being
    computed: ``singlet'', ``triplet'' or ``both''. This tag is only applicable
    for spin restricted calculation. For calculations in the ``triplet'' or
    ``both'' cases, \kw{SpinConstants} must be supplied (see
    p.~\pref{sec:dftbp.SpinConstants}).

  \item[\is{EnergyWindow}]\modif{\modtype{energy}} Energy range above the last
    transition at \is{NrOfExcitations} to be included in excited state spectrum
    calculation.

  \item[\is{OscillatorWindow}]\modif{\modtype{Dipole moment}} Screening cut-off
    below which single particle transitions are neglected in excitation spectra
    calculations. This selects from states above the top of the
    \is{EnergyWindow} (if present). This keyword should not be used if
    calculating forces or other excited state properties.

  \item[\is{WriteTransitions}] If set to \is{Yes}, the file TRA.DAT is
    created. This file contains a description of each requested excited state in
    terms of its single-particle transitions.

  \item[\is{WriteSPTransitions}] If set to \is{Yes}, the file SPX.DAT is
    created, which contains the spectrum at the uncoupled DFTB level (i.e.\ the
    single-particle excitations).

  \item[\is{WriteTransitionCharges}] For range separated calculations
    (p.~\pref{sec:dftbp.RangeSep}) setting this to \is{Yes} generates the file
    ATQ.DAT containing the Mulliken charges for a specific transition (chosen by
    \is{StateOfInterest}).

  \item[\is{WriteMulliken}] If set to \is{Yes}, the files XCH.DAT and XREST.DAT
    are created. The former contains atom-resolved Mulliken (gross) charges for
    the excited state of interest, the latter the excited-state dipole moment of
    the state.

  \item[\is{WriteCoefficients}] If set to \is{Yes}, the file COEF.DAT is
    created. This file contains the complex eigenvectors (molecular orbital
    coefficient) for the excited state of interest. They are derived from the
    relaxed excited state density matrix.

  \item[\is{WriteEigenvectors}] If set to \is{Yes}, the file excitedOrbs.bin is
    created (and excitedOrbs.out depending on whether \is{EigenvectorsAsText} is
    true). This file contains the natural orbitals for the specified excited
    state.

  \item[\is{WriteDensityMatrix}] If set to \is{Yes}, the file(s) DM?.bin are
    created. These hold the excited state density matrices, in the atomic
    orbital basis, for the evaluated states.

  \item[\is{TotalStateCoeffs}] Option to control data from
    \is{WriteCoefficients} or \is{WriteEigenvectors}. If set to \is{No} the
    total charge density of the output orbitals corresponds to the change in
    charge from the ground to excited state. If set to \is{Yes} instead it
    corresponds to the total charge density in the excited state.

  \item[\is{WriteXplusY}] If set to \is{Yes}, the file XplusY.DAT is
    created. This file contains the RPA vector $(X+Y)^{I\Sigma}_{ia}$ for all
    excited states (c.f., Eqn.~(18) in Ref.~\cite{heringer2007aes}).

  \item[\is{WriteTransitionDipole}] If set to \is{Yes}, the file TDP.DAT is
    created. This file contains the Mulliken transition dipole for each excited
    state.

  \item[\is{WriteStatusArnoldi}] If set to \is{Yes}, the file ARPACK.DAT is
    created, which allows the user to follow the progress of the Arnoldi
    diagonalisation.

  \item[\is{TestArnoldi}] If set to \is{Yes}, the file TEST\_ARPACK.DAT is
    created, which gives data on the quality of the resulting eigenstates.

  \item[\is{ExcitedStateForces}] If set to \is{Yes}, evaluated forces include
    the contributions from an excited state of interest. By default, it is set
    to \is{Yes} if forces are being calculated (for example in geometry
    optimisation) and to \is{No} otherwise. By setting it explicitly to \is{No},
    you can calculate the excitations during a molecular dynamics simulation
    that is being driven by the ground state forces only.

  \item[\is{CacheCharges}] If set to \is{No}, transition charges are calculated
    on the fly during the excited states calculation, instead of being
    cached. This makes the calculation considerably slower, but can help to
    decrease memory use substantially, if you are short on memory.

  \end{description}

\subsection{PP-RPA}
\label{sec:dftbp.pprpa}

This tag contains the specifications for the calculation of excitation energies
using the particle-particle random phase approximation
(pp-RPA)~\cite{Yang2017}. This approach, unlike time-dependent DFTB, allows the
computation of double and charge-transfer transitions. However it has the
limitation that the computed excitations have to involve, at least partially,
the highest occupied molecular orbital (HOMO) of the system.

For the computation of the excitation energies of a neutral $N$-electron system,
one needs to set up a ground state calculation for the two-electron deficient
($N-2$) system, i.e.\ a net \mbox{\kw{Charge} = +2.0} calculation. Please note
that if \kw{Charge} is set to 0.0 (the default value), the obtained transition
energies will correspond to a net negative $-2$ charged system (i.e.\ $N+2$
electrons). The system of interest must be closed-shell, therefore the
calculation must also be spin-restricted and performed for an even number of
electrons.

The pp-RPA method is currently implemented only for SCC-DFTB level excitations,
but can be performed on top of both SCC-DFTB or range-separated reference ground
state calculations. The SCC-DFTB calculations can be performed for molecular or
gamma point periodic systems, but the range-separated calculations can only use
molecular boundary conditions.

The specifications for this block have the following properties:

\begin{ptable}
  \kw{NrOfExcitations}   & i & & -  & \\
  \kw{Symmetry}          & s & & -  & \\
  \kw{NrOfVirtualStates} & i & & 0  & \\
  \kw{TammDancoff}       & l & & No & \\
  \kw{HHubbard}          & p & & -  & \\
\end{ptable}

\begin{description}

\item[\is{NrOfExcitations}] Specifies the number of vertical excitation energies
  to be computed for each symmetry (singlet or triplet).

\item[\is{Symmetry}] Specifies the spin symmetry of the excited states being
  computed: \is{singlet}, \is{triplet} or \is{both}. Please note that,
  \is{triplet} and \is{both} are effectively similar, as both singlet and
  triplet excitation energies will be printed out if either of these keywords
  are used.

\item[\is{NrOfVirtualStates}] Optional orbital constraint to speed up the
  calculation. It specifies the number of virtual states entering the pp-RPA
  equation. If set to zero or greater than the total number of virtual states of
  the system no constraint will be applied.

\item[\is{TammDancoff}] If set to \is{Yes}, the Tamm-Dancoff approximation will
  be employed. This will speed up the calculation.

\item[\is{HHubbard}] Hubbard-like parameters for each atom type including only
  the Hartree kernel. Values of some of these parameters are given in
  appendix~\ref{app:hhubbard}.

\end{description}

The output of the pp-RPA calculation are described in
section~\ref{sec:ppRPAout}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Time propagation of electronic state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ElectronDynamics}

\label{sec:dftbp.ElectronDynamics}

Calculate the real-time propagation of both the electronic state of the system
and the nuclei within the Ehfenfest approximation in the presence of an optional
external time-dependent perturbation, as described in Ref.~\cite{realtime1}.
\begin{ptable}
  \kw{Steps} &i &  & - & \\
  \kw{TimeStep} & r & & - & \\
  \kw{FieldStrength} & r & & - & \\
  \kw{Perturbation} & m & & None & \pref{sec:dftbp.tdpert} \\
  \kw{EnvelopeShape} & m & Perturbation = Laser\cb & Constant & \pref{sec:dftbp.envelope} \\
  \kw{Populations} & l & & No & \\
  \kw{Restart} & l & & No & \\
  \kw{RestartFromAscii} & l & Restart = Yes & No & \\
  \kw{WriteRestart} & l & & Yes & \\
  \kw{WriteAsciiRestart} & l & WriteRestart = Yes & No & \\
  \kw{WriteFrequency} & i & & 50 & \\
  \kw{WriteEnergyAndCharges} & l & & & \pref{sec:dftbp.tdextra} \\
  \kw{RestartFrequency} & i & & Steps/10 & \\
  \kw{Forces} & l & & No & \\
  \kw{IonDynamics} & l & & No & \\
  \kw{InitialTemperature} & r & IonDynamics = Yes & - & \\
  \kw{MovedAtoms} & (i|s)+ & IonDynamics = Yes & 1:-1 & \\
  \kw{Velocities} & (3r)* & IonDynamics = Yes & - & \\
  \kw{Pump} & l & & No & \\
  \kw{PumpProbeFrames} & i & Pump = Yes & - & \\
  \kw{ProbeProbeRange} & 2r & Pump = Yes & 0 Steps*TimeStep & \\
  \kw{Probe} & l & & No & \\
  \kw{EulerFrequency} & i & & 0 & \\
  \kw{WriteBondEnergy} & l & & No & \\
  \kw{WriteBondPopulation} & l & & No & \\
  \kw{FillingsFromFile} & l & & No & \\
\end{ptable}

\begin{description}
\item[\is{Steps}] Number of propagation steps to perform.

\item[\is{TimeStep}]\modif{\modtype{time}} Time interval between two
  electronic propagation steps.

\item[\is{FieldStrength}]\modif{\modtype{Electric field strength}}
  Peak intensity of the applied time-dependent electric field.

\item[\is{Populations}] If time-dependent (ground state) molecular
  orbital populations should be calculated. This is done by projecting
  the electronic state onto the ground state molecular orbitals. If
  Yes, they are saved every \is{WriteFrequency} steps to
  \verb|molpopul1.dat| (and \verb|molpopul2.dat| is the calculation is
  collinearly spin polarised).

\item[\is{Restart}] If calculation should be restarted. A restart file
  \verb|tddump.bin| (or \verb|tddump.dat|) must be available in the working
  directory.

\item[\is{RestartFromAscii}] If this is a restart calculation (\is{Restart} =
  \is{Yes}), read \verb|tddump.dat| instead of \verb|tddump.bin| in the working
  directory.

\item[\is{WriteRestart}] If a restart tddump file should be written. Details of
  the contents of the file are discussed in appendix~\ref{app:tddump}.

\item[\is{WriteAsciiRestart}] If a restart file is being written, select between
  producing \verb|tddump.bin| (\is{No}) or if set to \is{Yes},
  \verb|tddump.dat|.

\item[\is{RestartFrequency}] Number of steps every which the system's
  state should be writen to restart file.

\item[\is{WriteFrequency}] Number of steps every which the atomic
  charges and molecular orbital populations (if \is{Populations} =
  Yes) are written to file.

\item[\is{Forces}] If forces on all atoms should be calculated and
  written to \verb|forcesvst.dat|.

\item[\is{IonDynamics}] If nuclei should be propagated in the dynamics
  using the Velocity Verlet algorithm. If set to Yes, then either the
  \is{InitialTemperature} keyworkd or the \is{Velocities} keyword must
  be set. The positions and velocities every \is{WriteFrequency} steps
  will be saved in XYZ format to \verb|tdcoords.xyz|.

\item[\is{MovedAtoms}] List of atoms that will be moved during the
  dynamics. See \is{MovedAtoms} in section
  \ref{sec:dftbp.SteepestDescent}.

\item[\is{InitialTemperature}]\modif{\modtype{energy}} Create starting
  velocities for the Ehrenfest MD according to the Maxwell-Boltzmann
  distribution at the specified temperature. This keyword is redundant in the
  case of specified initial velocities or for a restarted trajectory, so is
  omitted in those cases.

\item[\is{Velocities}]\modif{\modtype{velocity}} Specified atomic
  velocities for all the atoms of the given structure (including
  ``velocities'' for any stationary atoms, which are silently
  ignored). (\textbf{Note:} if the velocities from a previous (MD or
  Ehrenfest) run are used, the velocities printed in the XYZ files are
  specified in {\AA}/ps, so this should be set in the input). See
  section \ref{sec:dftbp.VelocityVerlet}.

\item[\is{Pump}] If this trajectory corresponds to the dynamics under
  a pump pulse, with the intention to probe the system afterwards at
  different delay times, to simulate a pump-probe transient absorption
  experiment. The effect of the keyword is to write dump files named
  \textit{i}\verb|ppdump.bin| ($i=0,\dots,\mathrm{PumpProbeFrames}$)
  containing the state of the system (density matrix, coordinates,
  velocities) every a given number of steps.

\item[\is{PumpProbeFrames}] Number of total snapshots of the system in
  the pump trajectory that will be saved for future probe
  simulations. These are spaced uniformly in time over the
  \is{PumpProbeRange}.

\item[\is{PumpProbeRange}]\modif{\modtype{time}} The time range
  (initial time and final time, separated by a space) between which
  the snapshots will be dumped to file.

\item[\is{Probe}] If this is the simulation of a probe. If set to yes,
  then automatically \is{Restart} = Yes, \is{WriteRestart} = No, and
  only the dipole moment output files will be written, since this
  keyword is used together with a \is{Kick} perturbation to probe the
  system. Notice that the dump file that contains the state of the
  system that you want to be probed must be present and renamed to
  \verb|tddump.bin| (or \verb|tddump.dat| if using the ascii
  version of this data).

\item[\is{EulerFrequency}] Number of steps every which an Euler
  integration step is done in the electronic propagation (that
  normally uses the Leapfrog algorithm). The default value ensures no
  Euler steps are done during the dynamics. If used, must be set to
  something larger than 50.

\item[\is{WriteBondEnergy}] If the pairwise non-self-consistent part of the bond
  energy, $BE(t)$, should be written to the \verb|bondenergy.bin| file every
  \is{WriteFrequency} steps. The populations are evaluated as
  \begin{equation*}
    BE_{AB}(t) = \sum_{\mu \in A} \sum_{\nu \in B} \rho_{\mu \nu}(t) H^0_{\mu
      \nu}(t),
  \end{equation*}
  where $A$ and $B$ are the atoms in the central cell and $S$ is the overlap
  matrix. They are stored in binary format (see \ref{sec:beformatout}).

\item[\is{WriteBondPopulation}] If the pairwise time-dependent bond
  Mulliken-like population, \index{bond populations} $BP(t)$, should be written
  to the \verb|bondpop.bin| file every \is{WriteFrequency} steps. The
  populations are evaluated as
  \begin{equation*}
    BP_{AB} (t) = \sum_{\mu \in A} \sum_{\nu \in B} \rho_{\mu \nu}(t) S_{\mu
      \nu}(t),
  \end{equation*}
  where $A$ and $B$ are the atoms in the central cell and $S$ is the overlap
  matrix. They are stored in binary format (see \ref{sec:bpformatout}).

\item[\is{FillingsFromFile}] If the initial fillings (molecular
  orbital occupations) should be read from a file named
  ``fillings.ini'' present in the working directory. If used, the
  ground state fillings are replaced by these and used to build the
  initial density matrix.

\end{description}

\subsection{\kwcb{Perturbation}}
\label{sec:dftbp.tdpert}

Determines the type of perturbation that is applied.

\paragraph{\iscb{None}}
No time-dependent perturbation is included, free dynamics are calculated.

\paragraph{\kwcb{Kick}}
Perform Dirac-delta perturbation (or a \textit{kick}) to the density matrix.

\begin{ptable}
  \kw{PolarizationDirection} & s &  & - & \\
  \kw{SpinType} & s & SpinPolarisation = Colinear \cb & Singlet & \\
\end{ptable}

\begin{description}
 \item[\is{PolarizationDirection}] The cartesian axis for the kick:
   \verb|"x"|, \verb|"y"| or \verb|"z"|. If set to \verb|"all"|,
   calculates the three directions $x,y,z$ consecutively. For
   polarization direction \verb|x,y,z| the dipole moment output file
   will be named \verb|mux.dat|, \verb|muy.dat|, \verb|muz.dat|,
   respectively.

 \item[\is{SpinType}] Must be eiher ``Singlet'' or ``Triplet'' for
   singlet or triplet spectra, respectively. Only implemented for
   collinear spin polarisation.
\end{description}

\paragraph{\kwcb{Laser}}
\label{sec:dftbp.tdpert.laser}
Apply time-dependent sinusoidal perturbation to the system.

\begin{ptable}
  \kw{PolarizationDirection} & 3r &  & - & \\
  \kw{ImagPolarizationDirection} & 3r & & 0 0 0 & \\
  \kw{LaserEnergy} & r & & & \\
  \kw{Phase} & r & & 0 & \\
  \kw{ExcitedAtoms} & (i|s)+ & & 1:-1 & \\
\end{ptable}

\begin{description}
 \item[\is{PolarizationDirection}] Vector along which the electric
   field is polarized.

 \item[\is{ImagPolarizationDirection}] Imaginary part of the
   polarization vector. Useful for circularly polarized fields.

 \item[\is{LaserEnergy }] \modif{\modtype{energy}} Energy $\hbar
   \omega$ of the laser.

 \item[\is{Phase}] Optional initial phase of laser field in radians,
   such that sinusoidal component of the field is described by
   $\sin(\omega t + \phi)$.

 \item[\is{ExcitedAtoms}] List of atoms that will be excited by the
   laser. The format is the same as the \is{MovedAtoms} keyword.

\end{description}

\textbf{Note:} when working with periodic systems, {\bf the
  polarization direction must be orthogonal to the periodic
  directions} of the system. Therefore, the code does not work with 3D
periodic systems with laser perturbations.

\paragraph{\kwcb{KickAndLaser}}
Apply a \is{Kick} plus a \is{Laser} to the system. Useful for probing
the excited state on the system while being driven by a laser. The
keywords are a combination of the ones for kick and laser.

\begin{ptable}
  \kw{KickPolDir} & i &  & - & \\
  \kw{LaserPolDir} & 3r &  & - & \\
  \kw{LaserImagPolDir} & 3r & & 0 0 0 & \\
  \kw{LaserEnergy} & r & & - & \\
  \kw{LaserStrength} & r & & - & \\
  \kw{Phase} & r & & 0 & \\
  \kw{ExcitedAtoms} & (i|s)+ & & 1:-1 & \\
\end{ptable}

\begin{description}
 \item[\is{KickPolDir}] Same as PolarizationDirection in \is{Kick}.

 \item[\is{SpinType}] Same as \is{Kick}.

 \item[\is{LaserPolDir}] Same as PolarizationDirection in \is{Laser}.

 \item[\is{LaserImagPolDir}] Same as ImagPolarizationDirection in \is{Laser}.

 \item[\is{LaserEnergy }] \modif{\modtype{energy}} Same as in \is{Laser}.

 \item[\is{LaserStrength}] Peak intensity of the applied laser (in
   this mode, \is{FieldStrength} is the kick intensity).

 \item[\is{Phase}] Same as in \is{Laser}.

 \item[\is{ExcitedAtoms}] Same as in \is{Laser}.
\end{description}

\subsection{WriteEnergyAndCharges}
\label{sec:dftbp.tdextra}

This enables the write out of additional data during the propagation. By
default, excitations using \is{Kick} or \is{KickAndLaser} excitations sets this
to false, while \is{Laser} or \is{None} generates the extra files. This option
controls write out of Mulliken charges during the calculation (\verb|qvst.dat|)
and energy components (\verb|energyvst.dat|). If forces are evaluated this
keyword also enables writing of (\verb|forcesvst.dat|) or if ion dynamics are
allowed the atom coordinates during the calculation (\verb|tdcoords.xyz|).

\subsection{\kwcb{EnvelopeShape}}
\label{sec:dftbp.envelope}

Determines the envelope $f(t)$ of the laser, such that the laser field
is $E(t) = E_0 f(t) \sin(\omega t + \phi) $

\paragraph{\kwcb{Constant}}

Constant envelope and equal to $f(t) = 1$. Produces continuous wave laser.

\paragraph{\kwcb{Gaussian}}

Applies a Gaussian envelope function $f(t) = \exp (-(t-t_m)^2 /
\beta^2)$, where $t_m$ is the time at which the pulse is centered and
$\beta = \tau / 2\sqrt{\pi}$, $\tau$ being the duration of the pulse.

\begin{ptable}
  \kw{Time0} & r &  & 0 & \\
  \kw{Time1} & r & & - & \\
\end{ptable}

\begin{description}
 \item[\is{Time0}] Time at which the pulse starts.
 \item[\is{Time1}] Time at which the pulse ends. Is equal to Time0 + $\tau$.
\end{description}


\paragraph{\kwcb{Sin2}}

Applies a sin$^2$ envelope function:
\begin{equation*}
  f(t) =
  \begin{cases}
    \sin^2 (\pi(t-t_0)/\tau) & t_0 \leq t \leq t_0+\tau \\ 0 & t < t_0 \text{ or
    } t > t_0 + \tau.\\
  \end{cases}
\end{equation*}

The properties of this method are the same as for the Gaussian
(\kw{Time0}, \kw{Time1}) and have the same meaning.




  Example:\invparskip
\begin{verbatim}
  ElectronDynamics = {
    Steps = 40000
    TimeStep = 0.1
    # Total time will be then 4000 a.u. = 96.8 fs
    FieldStrength [v/a] = 0.01
    Perturbation = Laser {
      PolarizationDirection = 0.5 0.5 0
      LaserEnergy [ev] = 2.55
    }
    EnvelopeShape = Sin2 {
      Time1 [fs] = 30.0
    }
    Populations = Yes
    IonDynamics = Yes
    InitialTemperature [k] = 0.0
    Pump = Yes
    PumpProbeFrames = 1000
    PumpProbeRange [fs] = 0.0 50.0
    EulerFrequency = 200
  }
\end{verbatim}

%\subsection{Perturbation\cb}
%\label{sec:dftbp.tdreaddm}
%
%Reads in a starting density matrix (or mixture of density matrices)
%file(s). These are compatible with files generated by \is{WriteDensityMatrix}.
%
%\begin{ptable}
%  \kw{Name} & s &  & "" & \\
%  \kw{Weight} & s & & 1.0 & \\
%\end{ptable}
%
%\begin{description}
% \item[\is{Name}] Name of binary file holding the density matrix.
% \item[\is{Weight}] If more than one file is specified, this is required to
%   weight the contributions from each density matrix.
%\end{description}
%
%
%For example, a mixture of the lowest two degenerate bright excited states of the
%\cf{C6H6} molecule:
%\begin{verbatim}
%ReadDensityMatrix = {
%    File = {
%        Name = "DM7.dat"
%        Weight = 0.5
%    }
%    File = {
%        Name = "DM8.dat"
%        Weight = 0.5
%    }
%}
%\end{verbatim}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  REKS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{REKS}
\label{sec:dftbp.REKS}

This block collects some options to calculate REKS in the context of DFTB.
The \kw{Reks} keyword expects either \is{None} (default -- no use of REKS calculation)
or the \is {SSR22}\cb\ block as value.

\begin{ptable}
  \kw{SSR22} & p & SCC = Yes, SpinPolarisation $=$ \{\}, SpinConstants $\neq$ \{\} & None & \\
\end{ptable}

\subsection{SSR22}
\label{sec:dftbp.SSR22}

This tag contains the specifications for a DFTB/SSR(2,2) calculation~\cite{Lee_JCTC_2019}, based on ensemble DFT theory.

\textbf{Note:} the \dftbp{} binary can be compiled with OpenMP (not MPI) parallelization
and DFTB/SSR calculation is not compatible with time-dependent DFTB calculation.
In addition, it is not compatible with spin-polarisation, but it requires spin constants to treat open-shell microstates.

In general, REKS calculation can be classified as single-state REKS, SA-REKS and SI-SA-REKS.
In single-state REKS, only ground state is calculated and it can treat the state with multireference
character. SA-REKS and SI-SA-REKS can calculate the vertical excitation energies. The difference
is that the state-interaction term is considered in SI-SA-REKS so that more accurate states can
be generated. The corresponding oscillator strengths as well as excited state geometry optimisation
can be performed with these options, details of the resulting output files are given in appendix
\ref{sec:reks_files}.

In the context of DFTB, current REKS calculation is compatible with following functionalities.
The range-separated functional, external point charges and dispersion corrections can be
calculated with single-state REKS, SA-REKS or SI-SA-REKS. For the periodic system, only
gamma point sampling is supported with REKS. Especially, the stress evaluation and lattice
optimisation is possible with only single-state REKS.
The specifications for this block have the following properties:

\begin{ptable}
  \kw{Energy} & m & & \cb & \\
  \kw{TargetState} & i & & 1 & \\
  \kw{TargetMicrostate} & i & \parbox{0.3\textwidth}{Functional $\neq$ \{ "PPS" \}, StateInteractions = No} & 0 & \\
  \kw{ReadEigenvectors} & l & & No \\
  \kw{FonMaxIter} & i & & 20 & \\
  \kw{Shift} & r & & 0.3 & \\
  \kw{SpinTuning} & $(\text{r})^*$ & & \cb & \\
  \kw{TransitionDipole} & l & \parbox{0.3\textwidth}{Functional $\neq$ \{ "PPS" \}, TargetMicrostate = 0} & No & \\
  \kw{Gradient} & m & & ConjugateGradient \cb & \\
  \kw{RelaxedDensity} & l & & No & \\
  \kw{NonAdiabaticCoupling} & l & \parbox{0.3\textwidth}{Functional $\neq$ \{ "PPS" \}, StateInteractions = Yes} & No & \\
  \kw{VerbosityLevel} & i & & 1 & \\
\end{ptable}

\begin{description}
\item[\is{Energy}] Choice of energy evaluation in REKS method. This \is{Energy} block has following options:

  \begin{ptable}
    \kw{Functional} & $(\text{s})^*$ & & & \\
    \kw{IncludeAllStates} & l & & No & \\
    \kw{StateInteractions} & l & & No & \\
  \end{ptable}

  \begin{description}
  \item[\is{Functional}] Specifies the minimized energy functional in DFTB/SSR.
    This keyword reads a block consisted of the energy functionals that you want to include in the calculation.
    In DFTB/SSR(2,2), there are two possible choices for the minimzied energy functionals.
    One is PPS and the other is (PPS+OSS)/2. The former represents single-state REKS and the latter
    shows SA-REKS or SI-SA-REKS. The detailed form of the block is shown in below examples.
    The inclusion of state-interaction terms is determined by \is{StateInteractions}.

  \item[\is{IncludeAllStates}] If set to \is{Yes}, all computable energy states from current energy functionals
    are included for SA-REKS or SI-SA-REKS calculations. When you calculate single-state REKS, this option
    does not affect the result of calculation. The PPS and OSS states are calculated when this option sets to \is{No},
    while the additional DES state can be included if this sets to \is{Yes}. If you want to the doubly-excited
    configuration, please set to \is{Yes}. The detailed explanation about PPS, OSS and DES states
    is given in the Ref.~\cite{Lee_JCTC_2019}

  \item[\is{StateInteractions}] If set to \is{Yes}, the state-interaction terms between SA-REKS states is included,
    thus it generates SI-SA-REKS states. In general, SI-SA-REKS state can provide more reliable state
    when you want to compute the excited states.
  \end{description}

\item[\is{TargetState}] Specifies the target state that should be calculated. These are numbered
  from the ground state as 1, and so on. Note that the ordering of this option is different with the option
  \is{StateOfInterest} in time-dependent DFTB calculation.

\item[\is{TargetMicrostate}] Specifies the target microstate that should be calculated. The electronic
  configuration is given in the Ref.~\cite{Lee_JCTC_2019} or the source code of \dftbp{}. In SSR(2,2),
  fifth microstate is triplet configuration, thus this microstate can be roughly considered as triplet state.

\item[\is{ReadEigenvectors}] If set to \is{Yes}, the initial molecular orbitals are read from the
  eigenvec.bin file. If not, the initial orbitals are obtained from the diagonalization of non-SCC Hamiltonian.

\item[\is{FonMaxIter}] Specifies the maximum number of iterations used in the optimization of
  fractional occupation numbers. In general, the value of 20 is enough to converge the fractional
  occupation numbers in SCC cycle.

\item[\is{Shift}] Specifies the level shift value used in SCC cycle. The shift value should be increased
  to converge the SCC cycle when the orbital energies of the active orbitals are close to each other.

\item[\is{SpinTuning}] Specifies the scaling constants for atomic spin constants. DFTB/SSR sometimes
  shows wrong spin contribution for triplet microstate, thus the scaling of atomic spin constants are needed
  to generate correct spin contribution for each microstate. The standard to determine the scaling constants
  is provided in the Ref.~\cite{Lee_JCTC_2019}. The number of elements of \is{SpinTuning} block becomes the number of
  atomic species, and the ordering of the elements is same as the ordering of atomic species in input geometry file.

\item[\is{TransitionDipole}] If set to \is{Yes}, the file tdp.dat is created. This file contains a
  description of transition dipole moment between the electronic states in SA-REKS or SI-SA-REKS.

\item[\is{Gradient}] Choice of gradient solver used in CP-REKS equations. This \is{Gradient} block has the following choices:

  \begin{description}
  \item[\is{ConjugateGradient}] Uses a congugate-gradient based algorithm. This
    algorithm is usually recommended since it is considerably faster than other algorithms.

    \begin{ptable}
      \kw{CGmaxIter} & i & & 20 & \\
      \kw{Tolerance} & r & & 1e-8 & \\
      \kw{Preconditioner} & l & & No & \\
      \kw{SaveMemory} & l & & No & \\
    \end{ptable}

    \begin{description}
    \item[\is{CGmaxIter}] Specifies the maximum number of iterations used in the conjugate-gradient
    based algorithm. In general, the value of 20 is enough to solve the CP-REKS equations.

    \item[\is{Tolerance}] Specifies the tolerance used in the conjugate-gradient based algorithm.

    \item[\is{Preconditioner}] If set to \is{Yes}, it uses a preconditiner in the conjugate-gradient based algorithm. In general,
      the convergence speed is increased when this option sets to \is{Yes}, thus this option is recommended.

    \item[\is{SaveMemory}] If set to \is{Yes}, some variables (an orbital hessian matrix and the H-xc kernel) which need large
      memory allocation are saved in the memory. If these variables are saved, then the computational
      speed also increases but it shows large memory allocation, increasing as $O(N_\mathrm{basis}^4)$.
      If set to \is{No}, the CP-REKS equations are solved without saving these variables, thus it is relatively
      slower than the case that you set to \is{Yes}. In general, \is{No} option is recommended for large systems.
    \end{description}

  \item[\is{Direct}] Uses a direct matrix-inversion multiplication algorithm. This algorithm is usually considerably slow.
  \end{description}

\item[\is{RelaxedDensity}] If set to \is{Yes}, the file relaxed\_charge.dat is created. This file contains
  a description of relaxed charges for \is{TargetState} or \is{TargetMicrostate}. The relaxed charges can
  be used with external point charges in QM/MM calculations.

\item[\is{NonAdiabaticCoupling}] If set to \is{Yes}, the nonadiabatic couplings between SI-SA-REKS
  states are calculated. This option cannot be used in single-state REKS or SA-REKS state.

\item[\is{VerbosityLevel}] Specifies the printing level in standard output. This option determines the output up to
  energy information (\is{VerbosityLevel} = 0), gradient information (\is{VerbosityLevel} = 1), detailed information about
  SCC cycle and timing in gradient calculation (\is{VerbosityLevel} = 2).
\end{description}

Example for 3state SI-SA-REKS calculation with nonadiabatic couplings and modified spin constants:
\begin{verbatim}
Reks = SSR22 {
  Energy = {
    Functional = { "PPS" "OSS" }
    IncludeAllStates = Yes
    StateInteractions = Yes
  }
  TargetState = 2
  TargetMicrostate = 0
  ReadEigenvectors = No
  FonMaxIter = 30
  Shift = 0.3
  SpinTuning = { 3.0 3.0 }
  Gradient = ConjugateGradient {
    CGmaxIter = 100
    Tolerance = 1.0E-8
    Preconditioner = Yes
    SaveMemory = Yes
  }
  RelaxedDensity = Yes
  NonAdiabaticCoupling = Yes
  VerbosityLevel = 1
}
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  ParserOptions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ParserOptions}
\label{sec:dftbp.ParserOptions}

This block contains the options, which are effecting only the
behaviour of the HSD parser and are not passed to the main
program.
\begin{ptable}
  \kw{ParserVersion} & i & & \textrm{current input version} & \\
  \kw{WriteHSDInput} & l & & Yes & \\
  \kw{IgnoreUnprocessedNodes} &l & & No & \\
  \kw{StopAfterParsing} &l & & No & \\
\end{ptable}
\begin{description}
\item[\is{ParserVersion}] Version number of the input parser, which the
  input file was written for. If you are using an input file, which
  was created for an older version of \dftbp{}, you should set it to
  the parser version number of that code version. (The parser version
  number is printed at the beginning of the program run to the
  standard output.) \dftbp{} internally converts the input to its
  current format. The processed input (written to \verb|dftb_pin.hsd|)
  is always in the current format, and the \is{ParserVersion} property
  in it is always set to be the current parser version.

  The parser version will be inferred in case \kw{InputVersion} is set
  in the level above, in this case \is{ParserVersion} cannot be specified.

\item[\is{WriteHSDInput}] Specifies, if the processed input should be
  written out in HSD format. (You shouldn't turn it off without
  really good reasons.)

\item[\is{IgnoreUnprocessedNodes}] By default the code stops if it
  detects unused or erroneous keywords in the input, which probably
  indicates error(s) in the input. This {\em dangerous} flag suspends
  these checks. Use only for debugging purposes.

\item[\is{StopAfterParsing}] If set to \is{Yes}, the parser stops
  after processing the input and written out the processed input to
  the disc. It can be used to make sanity checks on the input without
  starting an actual calculation.

\end{description}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Parallel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parallel}
\label{sec:dftbp.Parallel}

This block contains the options, which are effecting the parallel behaviour of
the code. They only take effect, if the code was compiled with MPI-support.

\begin{ptable}
  \kw{Groups} & i & 1 & & \\
  \kw{UseOmpThreads} & l & .false. & & \\
  \kw{Blacs} & p & \cb & & \\
\end{ptable}
\begin{description}
\item[\is{Groups}] Number of process groups. Specifying more than one process
  group enables parallelisation over k-points and spin, as processes in
  different process groups are working on different k-points and spins at the
  same time. The number of process groups must be a divisor of the total number
  of MPI-processes. Default: 1 (all processes work at the same k-point and spin
  at a given time). Note that transport calculations between contacts are
  currently incompatible with multiple process groups (see
  section~\ref{app:transp}).

\item[\is{UseOmpThreads}] Enables the usage of OpenMP-threads (hybrid
  MPI/OpenMP-parallelisation). In order to prevent you from accidently running
  more processes and threads than appropriate for your hardware, this feature is
  turned off by default. Consequently in this case the MPI-parallelised binary
  will stop if the maximal number of OpenMP-threads is greater than one when
  DFTB+ is started. (You can usually set the number of maximally allowed
  OpenMP-threads by setting the \is{OMP\_NUM\_THREADS} environment variable in
  your shell.)

  You can enable this option if you wish to run DFTB+ with hybrid
  parallelisation. You would then typically start fewer MPI-processes than
  physical cores on each node and also set the number of threads accordingly.
  This is currently an experimental feature in DFTB+ and is recommended for
  experienced users only.

\item[\is{Blacs}] Contain BLACS specific settings. Currently only supports
  \is{BlockSize}, which specifies the row and column block size for the
  block-cyclic distributions (with default size of 32).

  Example:
\begin{verbatim}
  Parallel {
    Groups = 2
    Blacs {
      BlockSize = 64
    }
  }
\end{verbatim}

\end{description}


\chapter{Output of \dftbp}
\label{sec:dftbp.output}

This chapter contains the description of some of the output files of
\dftbp{} where the output format is not self documenting. Unless
indicated otherwise, numbers in the output files are given in atomic
units (with Hartree as the energy unit).

\section{band.out}
\label{sec:dftbp.bandout}
\index{band.out}

This contains the band energies and occupation of levels in electron volts and
electron charge units as columns one and two. The file is printed if
\is{WriteBandOut} = \is{Yes} (see section \ref{sec:dftbp.Analysis}). Blocks of
numerical results start with a line which labels the k-point and spin channel
for the energies.

See the \dptools{} package for utilities for converting the data in this file
into band-structures and density of states information suitable for plotting.

\section{detailed.out}
\label{sec:dftbp.detailedout}
\index{detailed.out}

This file contains details of the total energy and its components, as well as
optional information on forces, atomic charges and other properties. It is
intended for quick viewing, while values given to more significant figures are
available in results.tag.

Some of the information available in the file will also depend on the method
being used in the calculation. For example, not all electronic solvers make the
ground state electronic entropy available, hence only the internal energy would
be quoted. Similarly, while the free energy of the system which when
differentiated by atomic coordinates or boundary conditions gives the forces or
stresses (printed as \textit{Force related energy}) this is not currently
available for some types of non-equilibrium transport calculations.

Some of the common energy results printed in this file are:\\

\begin{center}
  \begin{tabular}{|l|p{5.5cm}|}
    \hline
    TS& Product of the electron entropy and temperature\\
    Total Electronic energy& The non-SCC energy plus other contributions to the
    electronic energy (SCC, spin, $\ldots$)\\
    Repulsive energy& The pairwise contribution to the total energy\\
    Total energy& Sum of electronic energy\\
    Extrapolated to 0& Estimated zero temperature energy if at finite
    temperatures\\
    Total Mermin free energy& $U - T S$, relevant free energy at finite
    temperatures\\
    Force related energy& Free energy relevant to forces in the system\\
    Gibbs free energy& Energy corrected by $- p V$, i.e. the pressure and volume\\
    MD Kinetic Energy& Kinetic energy of atoms in molecular dynamics\\
    Total MD Energy& Sum of finite temperature electronic, repulsive and atomic
    kinetic energies\\
    \hline
  \end{tabular}
\end{center}

Where available the Fermi level $\mu$ (i.e. the chemical potential of the
electrons in the system) is also printed. For systems with an externally fixed
Fermi level (i.e. where the total charge can change), this contribution is
included in the Force related energy:
\begin{equation*}
  \Delta E = + q_\mathrm{total} \mu,
\end{equation*}
but for calculations with fixed numbers of electrons it is not included in this
energy. \textbf{Note:} The total energy reference may not match your required
case in some situations, for example a shift with respect to the average
electrostatic potential (in periodic cases) or whether the chemical potential
should be with respect to the valence band maximum may be needed (see for
example the discussion in Ref.~\cite{Lany_2009}).

\section{results.tag}
\label{sec:dftbp.results}
\index{results.tag}

This contains machine readable results labeled with the type and size of the
data blocks. The results are given in atomic units and are formatted as:
\begin{verbatim}
label               :type:shape:
\end{verbatim}

The variable type is real, complex, integer or logical. The shape information is
\newline :ndim: size$_1$,size$_2$,$\ldots$,size$_{ndim}$:\newline where ndim is
the number of dimensions, organised with the Fortran convention and of size
size$_1$ $\times$size$_2$ $\times$size$_2$ $\times \ldots$.

In the special case of scalar variables the shape is :0:.

A typical example of mixed scalar and both one and two dimensional results would
be similar to:
\begin{verbatim}
mermin_energy       :real:0:
 -0.672967201447815E+000
total_energy        :real:0:
 -0.672879398682698E+000
forces              :real:2:3,3
 -0.243590222274811E+000 -0.199780753617099E-001 -0.000000000000000E+000
  0.465478448963764E+000 -0.228550455811745E+000 -0.000000000000000E+000
 -0.221888226688953E+000  0.248528531173455E+000 -0.000000000000000E+000
gross_atomic_charges:real:1:3
  0.171448741143825E+000 -0.254714832621691E+000  0.832660914778645E-001
\end{verbatim}

\section{hamsqrN.dat, oversqr.dat}
\label{sec:dftbp.hamsqr}
\index{hamsqr.dat}\index{oversqr.dat}

The files \verb|hamsqrN.dat| and \verb|oversqr.dat| contain the square
(folded) Hamiltonian and overlap matrices. The number \verb|N| in the
filename \verb|hamrealN.dat| indicates the spin channel. For spin
unpolarised calculation it is 1, for spin polarised calculation it is
1 and 2 for spin-up and spin-down, respectively while for non-collinear
spin it is charge, $x$, $y$ and $z$ for 1, 2, 3 and 4. Spin orbit is
not currently supported for this option.

Only non-comment lines (lines not starting with "\#") are documented:
\begin{itemize}

\item Flag for signalling if matrix is real (\verb|REAL|), number of
  orbitals in the system (\verb|NALLORB|), number of kpoints
  (\verb|NKPOINT|). For non-periodic (cluster) calculations, the
  number of kpoints is set to 1.

\item For every $k$-point:
  \begin{itemize}
  \item Number of the $k$-point. For molecular (non-periodic)
    calculations only 1 $k$-point is printed.
  \item The folded matrix for the given $k$-point. It consists of
    \verb|NALLORB| lines $\times$ \verb|NALLORB| columns. If the
    matrix is not complex (\verb|REAL| is \verb|F|), every column
    contains two numbers (real and imaginary part).
  \end{itemize}
\end{itemize}

The files are produced if requested by \is{WriteHS} = \is{Yes} (see
section~\ref{kw:dftbp.WriteHS}).

\section{hamrealN.dat, overreal.dat}
\label{sec:hamreal}
\index{hamreal.dat}\index{overreal.dat} The files \verb|hamrealN.dat|
and \verb|overreal.dat| contain the real space Hamiltonian and overlap
matrices. The number \verb|N| in the filename \verb|hamrealN.dat|
indicates the spin channel. For spin unpolarised calculation it is 1,
for spin polarised calculation it is 1 and 2 for spin-up and
spin-down, respectively, while for non-collinear spin it is charge,
$x$, $y$ and $z$ for 1, 2, 3 and 4. Spin orbit is not currently
supported for this option.

Note: The sparse format contains only the "lower triangle" of the real
space matrix. For more details about the format and how to obtain the
upper triangle elements, see reference~\cite{dftbp-2007paper}. Also note,
that for periodic systems the sparse format is based on the
\emph{folded} coordinates of the atoms, resulting in translation
vectors (ICELL) which look surprising at first glance.

Only non-comment lines (lines not starting with "\#") are documented:
\begin{itemize}
\item Number of atoms in the system (\verb|NATOM|)
\item For every atom:
  \begin{itemize}
  \item Atom number (\verb|IATOM|), number of neighbours including the
    atom itself (\verb|NNEIGH|), number of orbitals on the atom
    (\verb|NORB|)
  \end{itemize}
\item For every neighbour of every atom:
  \begin{itemize}
  \item Atom number (\verb|IATOM1|), neighbour number (\verb|INEIGH|),
    corresponding image atom to the neighbour in the central cell
    (\verb|IATOM2F|), coefficients of the translation vector between
    the neighbour and its corresponding image (\verb|ICELL(1)|,
    \verb|ICELL(2)|, \verb|ICELL(3)|). Between the coordinates of the
    neighbour $\mathbf{r}_{\text{INEIGH}}$ and the image atom
    $\mathbf{r}_{\text{IATOM2F}}$ the relation
    \begin{equation*}
      \mathbf{r}_{\text{INEIGH}} = \mathbf{r}_{\text{IATOM2F}} + \sum_{i=1}^3
      \text{ICELL}(i)\, \mathbf{a}_i
    \end{equation*}
    holds, where $\mathbf{a}_i$ are the lattice vectors of the supercell.
  \item The corresponding part of the sparse matrix. The data block
    consists of \verb|NORB(IAT1)| lines and \verb|NORB(IAT2F)| columns.
  \end{itemize}
\end{itemize}

The files are produced if requested by \is{WriteRealHS} = \is{Yes}
(see section~\ref{kw:dftbp.WriteRealHS}).

\section{eigenvec.out, eigenvec.bin}
\label{sec:dftbp.eigenvec}
\index{eigenvec.out}\index{eigenvec.bin}

These files contain the eigenvectors from the Hamiltonian, stored
either as plain text (eigenvec.out) or in the native binary format of
your system (eigenvec.bin).

The plain text format file \verb|eigenvec.out| contains a list of the values of
the components of each eigenvector for the basis functions of each atom. The
atom number in the geometry, its chemical type and the particular basis function
are listed, followed by the relevant value from the current eigenvector and then
the Mulliken population for that basis function for that level\index{state
  resolved Mulliken population}. The particular eigenvector, $k$-point and spin
channel are listed at the start of each set of eigenvector data. In the case of
non-collinear spin, the format is generalised for spinor wavefunctions. Complex
coefficients for both the up and down parts of the spinors are given (instead of
single eigenvector coefficient) followed by four values -- total charge, then
$(x,y,z)$ magnetisation.

The binary format file \verb|eigenvec.bin| contains the (unique) runId of the
DFTB+ simulation which produced the output followed by the values of the
eigenvectors. The eigenvector data is ordered so that the individual components
of the current eigenvector are stored, with subsequent eigenvectors for that
$k$-point following sequentially. All $k$-points for the current spin channel
are printed in this order, followed by the data for a second channel if spin
polarised.

The files are produced if requested by setting \is{WriteEigenvectors} =
\is{Yes}, with \is{EigenvectorsAsText} being also required to produce the plain
text file (see section~\ref{kw:dftbp.WriteEigenvectors} for details).

\section{charges.bin / charges.dat}
\label{sec:charges.bin}
\index{charges.bin}

The file \verb|charges.bin| contains the orbitally-resolved charges for each
atom. In later versions of \dftbp{} this format includes a check sum for the
total charge and magnetisation. In the case of orbital potentials
(p.~\pref{sec:DFTB+U}) the file also contains extra population information for
the occupation matrices.

This file is produced as part of the mechanism to restart ground state SCC
calculations, see sections~\ref{kw:dftbp.RestartFrequency}
and~\ref{kw:dftbp.MDRestartFrequency}.

Equivalent data can also be present in the file \verb|charges.dat|, but stored
as plain text. The options \is{WriteChargesAsText} and \is{ReadChargesAsText}
control which cases are generated and read respectively.

Appendix \ref{app:restartfiles} contains details of the contents of the file.

\section{md.out}
\label{sec:md.out}
\index{md.out}

This file is only produced for \iscb{VelocityVerlet} calculations (See
p.~\pref{sec:dftbp.VelocityVerlet}). It contains a log of information generated during MD
calculations, and appended every \kw{MDRestartFrequency} steps. In the case of
small numbers of atoms and long MD simulations it may be useful to set
\is{WriteDetailedOut} to \is{No} and examine the information stored in this file
instead.

\section{Electrostatic potential data}
\label{sec:dftbp.espfile}
\index{ESP.dat}

The output from evaluating the electrostatic potential (see
page~\pref{sec:dftbp.esp}). The first line consists of a comment mark followed
by a logical variable as to whether there is an external electric field (or
not), followed by 3 values for any regular grid pattern present in the system
and the total number of points. If the data is gridded, the next four lines
contain the origin and grid separation vectors in {\AA}ngstroms.

The next line is a comment, then the locations and the potential experience for
a positive charge due to the internal field plus optionally the external field
(from point charges or homogeneous electric fields). Values are given in
Volts. In the case of gridded data, the location field is omitted.

For an example with a regular grid
\begin{verbatim}
#  T     1     1     1 1
#  0.000000000000E+00 -0.200000000000E+01 -0.200000000000E+01
#  0.200000000000E+01  0.000000000000E+00  0.000000000000E+00
#  0.000000000000E+00  0.200000000000E+01  0.000000000000E+00
#  0.000000000000E+00  0.000000000000E+00  0.200000000000E+01
# Internal (V)        External (V)
  0.173386318927E-10  0.314737193575E+00
\end{verbatim}

In the case where there is no regular grid:
\begin{verbatim}
#  T     0     0     0 1
#           Location (AA)             Internal (V)        External (V)
  0.0000E+00 -0.2000E+01 -0.2000E+01  0.173386318927E-10  0.314737193575E+00
\end{verbatim}

In the case where data is generated for multiple geometry steps, this is also
shown in the label:
\begin{verbatim}
#  F     1     5     5 25
#  0.000000000000E+00 -0.200000000000E+01 -0.200000000000E+01
#  0.100000000000E+01  0.000000000000E+00  0.000000000000E+00
#  0.000000000000E+00  0.100000000000E+01  0.000000000000E+00
#  0.000000000000E+00  0.000000000000E+00  0.100000000000E+01
# Internal (V) Geo 0
  0.215249473376E-01
  .
  .
# Internal (V) Geo 10
  0.215815549672E-01
  .
  .
\end{verbatim}

\section{Excited state results files}
\label{sec:tddftb_lr}
Several files are produced during excited state calculations depending on the
particular settings from section~\ref{sec:dftbp.ExcitedState}.

\textbf{Note:} in the case of degeneracies, the oscillator strengths depend on
arbitrary phase choices made by the ground state eigensolver. Only the sum over
the degenerate contributions is well defined for most single particle transition
properties, and label ordering of states may change if changing eigensolver or
platform. For the excited state, properties like the intensities for
individual excitations in degenerate manifolds again depend on phase choices
made by both the ground and excited eigensolvers.


\subsection{ARPACK.DAT}
\index{ARPACK.DAT}

Internal details of the ARPACK solution vectors, see the ARPACK
documentation~\cite{Lehoucq97arpackusers} for details.

\subsection{COEF.DAT}
\index{COEF.DAT}

Data on the projection of this specific excited state onto the ground state
orbitals. For the specific exited state, the (complex) decomposition of its
single particle states onto the ground state single particle levels, together
with its fractional contribution to the full excited state are given.

General format:

\begin{tabular}{|l|p{5.5cm}|}
\hline
{\small T F                                }&Legacy flags\\
{\small   1  1.9999926523  2.0000000000    }&level 1, fraction of total WF, 2.0\\
{\small-0.1944475716  0.0000000000 -0.1196876988  0.0000000000 ....    }&real then
imaginary projection of level 1\\
                                    &onto ground state 1, then ground state 2, etc.\\
{\small-0.1196876988  0.0000000000 -0.1944475703  0.0000000000 ....    }&\\
{\small.}&\\
{\small.}&\\
{\small.}&\\
{\small   2  1.9999866161  2.0000000000    }&level 2\\
{\small-0.2400145188  0.0000000000 -0.1767827333  0.0000000000 ....}& real then
imaginary projection of state 2\\
{\small.}&\\
{\small.}&\\
{\small.}&\\ \hline
\end{tabular}

\subsection{EXC.DAT}
\index{EXC.DAT}

Excitations data including the energies, oscillator strength, dominant Kohn-Sham
transitions and the symmetry.

Example first few transitions for C$_4$H$_4$:

\begin{verbatim}
     w [eV]       Osc.Str.         Transition         Weight      KS [eV]    Sym.

 =========================================

      5.551       0.5143882       11   ->    12        1.000       4.207      S
      5.592       0.0000000       10   ->    12        1.000       5.592      S
\end{verbatim}

Two examples of singlet transitions with energies of 5.551 and 5.592~eV. The
first is dipole allowed, the second not. In both cases they are transitions
primarily (weight of 1.000) to single particle state 12, and are of singlet
character (``S'').

In the case of spin-polarised calculations, an additional column of values are
given instead of the symmetry, showing the level of spin contamination in the
state (labelled as \verb|D<S*S>|), with typically states where a magnitude of
less than 0.5 is usually considered reliable~\cite{garcia14Thesis}.

\subsection{SPX.DAT}
\index{SPX.DAT}

Single particle excitations (SPX) for transitions between filled and empty single
particle states of the ground state. These are given in increasing single
particle energy and show the oscillator strength and index of the Kohn-Sham-like
states that are involved.

\begin{verbatim}
       #      w [eV]       Osc.Str.        Transition

 ============================

       1      5.403       0.2337689       15   ->    16
       2      5.403       0.2337689       14   ->    16
       3      5.403       0.2337689       15   ->    17
       4      5.403       0.2337689       14   ->    17
       5      6.531       0.0000000       13   ->    16
       6      6.531       0.0000000       12   ->    16
\end{verbatim}

\subsection{TDP.DAT}
\index{TDP.DAT}

Detail of the magnitude and direction of the transition dipole from the ground
to excited states.

\subsection{TRA.DAT}
\index{TRA.DAT}

Decomposition of the transition from the ground state to the excited states. The
energy and spin symmetry are given together with the contributions from each of
the single particle transitions.

\subsection{TEST\_ARPACK.DAT}
\index{TEST\_ARPACK.DAT}

Tests on the quality of the eigenvalues and vectors returned by ARPACK. For the
$i^\mathrm{th}$ eigen-pair, the eigenvalue deviation corresponds to the
deviation from $\left( \langle \mathbf{x}_i | H | \mathbf{x}_i\rangle -
\epsilon_i \right)$, The eigen-vector deviation is a measure of rotation of the
vector under the action of the matrix: $\left| \left( H | \mathbf{x}_i\rangle -
\epsilon_i | \mathbf{x}_i\rangle \right) \right|_2$, the normalisation deviation
is $\langle \mathbf{x}_i | \mathbf{x}_i\rangle - 1$ and finally largest failure
in orthogonality to other eigenvectors is given.

Example:\\
\begin{tabular}{lllll}
{\tt State} & {\tt Ei deviation} & {\tt Evec deviation} & {\tt Norm deviation} &
{\tt Max non-orthog}\\ {\tt 1} & {\tt -0.19428903E-15} & {\tt 0.80601119E-15} &
{\tt 0.19984014E-14} & {\tt 0.95562226E-15}\\ {\tt 2} & {\tt 0.27755576E-16} &
{\tt 0.85748374E-15} & {\tt 0.48849813E-14} & {\tt 0.36924443E-15}\\ {\tt 3} &
{\tt -0.12490009E-15} & {\tt 0.88607302E-15} & {\tt 0.88817842E-15} & {\tt
  0.60384195E-15}\\
\end{tabular}

\subsection{XCH.DAT}
\index{XCH.DAT}

Charges on atoms in the specified excited state. The top line contains the
symmetry (Singlet or Triplet) and the number of the excited state. The next line
is the number of atoms in the structure followed by some header text. Then on
subsequent lines the number of each atom in the structure and its charge are
printed.

\subsection{XplusY.DAT}
\index{XplusY.DAT}

Expert file with the RPA  $(X+Y)^{I\Sigma}_{ia}$ data for all the calculated
excited states.

Line 1: number of single particle excitations and the number of calculated
excited states\\
Line 2: Level number 1, nature of the state (S, T, U or D) then excitation
energy (in Hartree)\\
Line 3: expansion in the KS single particle transitions\\
.\\
.\\
.\\
Line 2: Level number 2, nature of the state (S, T, U or D) then excitation
energy (in Hartree)\\

\subsection{XREST.DAT}
\index{XREST.DAT}

Dipole moment of the specified excited state in units of Debye.

\section{Electron dynamics results files}
\label{sec:tddftb_rt}

Real-time dynamics simulations produce the following output files:

\subsection{energyvst.dat}
\index{energyvst.dat}

Time-dependent energy components in Hartree. The column order is the
following:

time (in fs), total energy, non-SCC energy, SCC energy, spin energy,
external field energy, repulsive energy, nuclear kinetic energy,
dispersion energy

\subsection{qsvst.dat}
\index{qsvst.dat}

Net atomic charges in electrons (negative atomic populations) as a
function of time, written every \kw{WriteFrequency} steps, with the
following column order:

time (in fs), net total charge, charge (atom 1), charge (atom 2), ...,
charge (atom N)

where N is the number of atoms, and the net total charge should always
be zero (up to numerical presicion) for a neutral system.

\subsection{mu.dat/mux.dat/muy.dat/muz.dat}
\index{mu.dat}

Dipole moment cartesian components $\mu_i$ as a function of time. The
name depends on the type or perturbation: for kicks the polarization
direction of the Dirac-delta field is indicated in the name of the
file; for lasers and pulses, the name is always mu.dat. The column
order for a spin unpolarised calculation is the following:

time (in fs), $\mu_x$, $\mu_y$, $\mu_z$

For (collinear) spin polarised calculations, the dipole components are
also spin-dependent:

Time (in fs), $\mu_x$ (up), $\mu_y$ (up), $\mu_z$ (up), $\mu_x$
(down), $\mu_y$ (down), $\mu_z$ (down).

\subsection{laser.dat}
\index{laser.dat}

Created if \kw{Perturbation} = \kwcb{Laser}
(see~\ref{sec:dftbp.tdpert.laser}). Laser field components $E_i$ (given in
V/\AA) as a function of time, with data in the following order:

time (in fs), $E_x$, $E_y$, $E_z$

\subsection{molpopul1.dat/molpopul2.dat}

Created if \kw{Populations} = Yes. The number $i$ in molpopul$i$.dat
indicates the spin channel for spin polarised calculations (if spin
unpolarised, only molpopul1.dat is written). It contains the
populations projected onto the ground state molecular orbitals (GSMO),
written every \kw{WriteFrequency} steps. The projection is done by
calculating first the change of basis matrix $\Lambda$ from the
coefficients matrix of the GSMO, $\Lambda = C^{-1}$, and then
projecting the time-dependent density matrix of the system in the
atomic orbital basis:
\begin{equation*}
  \rho^{GSMO}(t) = \Lambda \rho^{AO}(t) \Lambda^{\dagger}
\end{equation*}
The diagonal elements of $\rho^{GSMO}$ are the sought populations. The first
column is the time (in fs), and from the second one the populations starting
from the lowest lying orbital.

\subsection{forcesvst.dat}

Time-dependent forces, calculated within the Ehrenfest approach and
ignoring velocity-dependent terms, printed every \kw{WriteFrequency}
steps, if \is{Forces} = Yes. The file has $3N +1$ columns where $N$ is
the number of atoms, in the following order:

time (in fs), $F_1^x$, $F_1^y$, $F_1^z$, $F_2^x$, $\dots$, $F_N^x$,
$F_N^y$, $F_N^z$

\subsection{tdcoords.xyz}

Coordinates and velocities of the atoms saved every
\kw{WriteFrequency} steps in XYZ format, if \is{IonDynamics} = Yes.

\subsection{bondenergy.bin}
\label{sec:beformatout}

Atom pair-wise non-self-consistent bond energy, written if \is{WriteBondEnergy}
= \is{Yes}, every \is{WriteFrequency} steps. All numbers are stored in the
following order (using the Fortran real precision that the code is compiled
with):

time (in fs), $BE_{1,1}$, $BE_{2,1}$, $BE_{3,1}$, $\dots$, $BE_{N,1}$,
$BE_{1,2}$, $BE_{2,2}$, $\dots$, $BE_{N,2}$, $\dots$, $BE_{N,N}$

\subsection{bondpop.bin}
\label{sec:bpformatout}

Atom pair-wise bond populations, written if \is{WriteBondOrder} = \is{Yes},
every \is{WriteFrequency} steps. (using the Fortran real precision that the code
is compiled with):

time (in fs), $BP_{1,1}$, $BP_{2,1}$, $BP_{3,1}$, $\dots$, $BP_{N,1}$,
$BP_{1,2}$, $BP_{2,2}$, $\dots$, $BP_{N,2}$, $\dots$, $BP_{N,N}$


\section{ppRPA\_ener.DAT}
\index{ppRPA\_ener.DAT}
\label{sec:ppRPAout}

Excitation energies obtained within the pp-RPA formalism (see
section~\ref{sec:dftbp.pprpa}). This output file also includes the
most dominant Kohn-Sham transition, its weight and energy difference
as well as the spin multiplicity of the excited state.

Here are, for instance, the first three singlet-singlet transitions
for furan:

\begin{verbatim}
     w [eV]        Transitions                         Weight            KS [eV]              Symm.

 =================================================

      6.411        HOMO -> LUMO + 0            0.998            5.162                     S
      6.904        HOMO -> LUMO + 1            0.973            6.755                     S
     11.339        HOMO -> LUMO + 0, 0        0.959            5.162,  5.162            S
\end{verbatim}

The first two excitations are single, whereas the third one is a
double transition with predominant HOMO-to-LUMO character.

\section{REKS results files}
\label{sec:reks_files}
Several files are produced during REKS calculations depending on the
particular settings from section~\ref{sec:dftbp.REKS}.


\subsection{tdp.dat}
\index{tdp-reks.dat}

Detail of the magnitude and direction of the transition dipole between all electronic states.

\subsection{relaxed\_charge.dat}
\index{relaxed\_charge.dat}

Charges on atoms in the specified state. The top line contains the total charge of the system.

%%% Local Variables:
%%% mode: latex
%%% TeX-master: "manual"
%%% End:
